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B2DE  - A Program  for  Solving  Systems  of  Nonlinear  Elliptic 
Partial  Differential  Equations  in  Two  Dimensions 

Manual  and  Users’  Guide 
James  L.  Blue 

Scientific  Computing  Division 

National  Bureau  of  Standards 
Gaithersburg,  MD  20899 

ABSTRACT 

B2DE  is  a program  for  solving  systems  of  nonlinear  elliptic  partial  differential 
equations  (PDEs)  in  two  dimensions.  The  program  is  a collection  of  modules  with 
an  interactive  driver.  Many  types  of  interactive  graphics  plots  are  included. 

Users  may  modify  the  driver,  and  may  be  able  to  construct  a “black  box” 
program  for  a restricted  class  of  PDEs. 

B2DE  is  available  from  the  author  at  the  above  address. 


1.  Introduction 

B2DE  is  a program  for  solving  systems  of  nonlinear  elliptic  partial  differential  equations 
(PDEs)  in  two  dimensions.  B2DE  has  been  used  for  solving  problems  as  simple  as  Laplace’s 
equation  [Kell86|  and  as  complicated  as  the  coupled  equations  of  electrostatics  and  chaxged- 
particle  transport  in  a transistor  [Blue83,  Wils85]. 

Since  there  is  no  algorithm  for  solving  all  solvable  nonlinear  systems  of  PDEs,  B2DE  cannot 
work  as  a “black  box”  program,  as  does,  for  example,  a Gaussian  elimination  program. 
Instead,  B2DE  is  a collection  of  modules;  the  user  must  decide  what  to  do  at  each  point 
of  the  solution  process.  An  interactive  driver  is  provided  for  the  user.  Expert  users  may 
wish  to  modify  the  driver,  and  may  construct  a “black  box”  program  for  a restricted  class 
of  PDEs.  Less  expert  users  may  wish  to  modify  some  routines  designed  for  easy  user 
modification. 

Several  kinds  of  graphic  displays  are  provided  to  aid  the  user  in  imderstanding  the  solution 
and  the  solution  process. 

This  manual  includes  a description  of  the  driver  program  and  a description  of  the  interface 
to  the  lower  level  programs.  It  does  not  include  a description  of  the  detailed  working  of 
the  lower  level  routines  themselves. 

The  drivers  and  example  routines  are  provided  in  well-structured  and  indented  Fortran  77. 
The  lower  level  routines  were  written  in  a private  version  of  Ratfor  [Kem75);  the  output 
of  the  Ratfor  processor  is  indented  Fortran  77,  but  is  not  pretty.  B2DE  is  single-precision. 
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Beginners  to  B2DE  should  first  try  one  of  the  sample  problems,  using  the  interactive 
driver.  Its  use  is  described  in  Section  3.2.  To  do  a new  problem,  users  must  modify  a 
main  program,  three  subroutines,  and  a data  file  to  describe  the  PDEs.  They  should  read 
Sections  3.1  - 3.3  and  Section  3.4.1. 

More  experienced  users  may  wish  to  change  some  of  the  USRxxx  routines.  They  should 
read  all  of  Section  3.4. 

Expert  risers  may  wish  to  modify  the  drivers,  the  DRVxxx  routines.  They  should  read 
Section  4. 

Installation  of  B2DE  on  your  computer  and  graphics  terminals  is  covered  in  Section  5. 
Section  6 contains  several  exzimple  problems. 

2.  What  B2DE  Does 
2.1  Mathematical  Model 

The  current  version  of  the  B2DE  solves  systems  of  N PDEs,  each  of  the  form: 


V-  [ai{x,y,u,du/dx,du/dy)'7ui]  = fi{x,y,u,duldx,du/dy)  i 


1,2,. ..,7V 


(1) 


in  a region  bounded  by  straight  line  segments  and  circular  arc  segments.  On  each  segment 
of  the  boimdary,  each  of  the  u*  must  obey  either  a nonlinear  Dirichlet  boundary  condition 
(type  1): 

ffi(2T,y,u)=0  (2) 

or  a nonlinear  normal  derivative,  of  Neumann,  boundary  condition  (type  2): 

dtL ' 

+ gi{x,y,u)  = 0.  (3) 

On  a given  boundary  segment,  the  type  of  boimdary  condition  may  be  different  for  dif- 
ferent PDEs.  This  is  a fairly  general  example  of  a system  of  elliptic  PDEs  arising  from 
conservation  of  a flux. 

The  numerical  software  was  inspired  by  finite-element  software  of  Bank  amd  Sherman 
[Bank79,  Bank82]  and  retains  most  of  their  philosophy. 

At  the  heart  of  the  package  is  a module  which  solves  a system  of  linearized  elliptic  PDEs 
using  linear  finite  elements  on  a mesh  of  triangles.  (Boundairy  “triangles”  may  have  one 
curved  side,  a circular  arc  segment.)  The  initial  triangulation  is  am  input  to  the  package; 
succeeding  triangulations  axe  calculated  adaptively  by  the  package.  For  a given  triangula- 
tion with  M vertices,  the  calculation  proceeds  as  follows. 

The  M finite-element  basis  functions  {6,^}  are  linear  on  eawrh  triangle;  bm  is  1 at  vertex 
m and  0 at  all  other  vertices.  For  the  tth  PDE,  the  set  of  vertices  on  type  1 boimdaries  is 
denoted  by  D{  (for  Dirichlet).  The  solution  is  approximated  by  a sum  of  basis  functions 

M 

ti.(a:,y)  = X]  “»m6m(x,y).  (•*) 

m=  1 
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The  coefficients  are  determined  by  a Galerkin  method  [Stra73];  the  error  in  the  ith  PDE 
is  made  orthogonal  to  eax:h  of  the  basis  functions  except  those  in  J9»-: 


j j [- V • (a,- Vu,)  + /,]  = 0,  m = 1, 2, . . . , M,  m not  in  Di  (5) 

The  remaining  conditions  on  the  a’s  are  that  the  type  1 boundary  conditions  hold  exactly: 

= 0,  m in  Di.  (6) 

Equation  (5)  is  integrated  by  parts  to  yield 

j j (a,  Vu,  • Vhm  + fihm)  + j cagibm.  = 0,  (7) 

where  the  single  integral  is  over  the  type-2  parts  of  the  boundary  only.  Equations  (6)  and 
(7)  comprise  MN  nonlinear  equations  in  MN  unknowns,  the  coefficients 

The  MN  nonlinear  equations  are  solved  by  an  iterative  process,  a damped  Newton’s 
method.  Let  v be  the  Newton  correction  to  u;  t;,-  has  the  same  form  as  u»-. 


M 


(8) 


n»=l 


To  calculate  v,  replace  u*  by  Ui  + in  (6)  and  (7)  and  linearize  in  v,-.  The  linearization 
of  (6)  yields 


N 


ym))  ^ ^ •_  j-»  fr\\ 

Qi  J/m))  ^ r Pjm  — Oj  771  in  2?,',  (9) 


and  the  linearization  of  (7)  yields 


jj{oiVvi 


^bm  + 


N 


Vu. . V6,„ 

y=i 


dai  dai  dvi  dai  dvi 


N r 


y=i 


duj  ^ d{dujldx)  dx  d{dujldy)  dx 
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dfi  dfi  dvj  dfi  dvj 

duj^^  d{dujldx)  dx  d{dujjdx)  dy 
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r ^ 

/‘■E 


dgi  dai  dai  dvj  dai  dvi 

+ Si—«i  + 9i -dv 0--  ; ^ + g.-  ' 


j=l  ^ 


duj  d{dujfdx)  dx  d{duj/dy)  dy 
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(10) 


Vu,-  • V6„»  + fib 


m 


Qi^m  — 


The  right-hand  side  is  the  residual,  r,m(u).  Substituting  for  u and  v from  (4)  and  (8) 
and  doing  the  integrals  numerically  give  a set  oi  MN  linear  equations  for  the  correction 
coefficients  {/3im}-  The  matrix  of  the  equations  is  the  Jacobian  matrix.  For  some  problems, 
the  Jstcobian  matrix  may  be  singular,  or  nearly  so.  B2DE  allows  the  user  to  “reguleLrize” 
the  Jacobian  matrix  by  adding  a constant  to  the  “a”  term  in  the  Jacobizin  matrix;  this 
is  effectively  like  adding  a multiple  of  the  identity  matrix  to  the  Jacobian  matrix.  (The 
residual  is  not  changed,  so  that  the  solution  foimd  is  the  correct  solution  to  the  nonlinear 
equations.) 

B2DE  allows  the  user  to  iterate  on  any  combination  of  the  PDEs,  leaving  the  others  fixed. 

There  is  no  guarantee  of  convergence  if  the  initial  guess  is  not  good  enough;  how  “good” 
is  “good  enough”  is  problem  dependent. 

These  linear  equations  are  sparse;  there  are  typically  about  7N  nonzero  elements  per  row, 
on  average,  and  MN  can  be  a few  himdred  or  a few  thousand.  On  the  initial  triangulation, 
with  MN  a few  hundred,  a satisfactory  way  to  solve  these  equations  is  directly,  with 
sparse  Gaussian  Elimination  [Sher78|.  On  later  triangulations,  with  MN  up  to  a few 
thouszmd,  direct  solutions  take  too  much  time  and  space;  iterative  methods  are  usually 
better,  especially  multi-level  iterative  methods. 

The  linear  equations  on  the  coarsest  level  are  solved  by  sparse  Gaussian  elimination.  The 
linear  equations  on  the  higher  levels  are  solved  by  multi-level  iteration  [Bank79,  Bank82, 
Bran77].  The  smoothing  iterations  provided  by  B2DE  are  Gauss-Seidel,  minimum-  residual 
Gauss-Seidel,  Block  Gauss-Seidel  using  N equations  per  block  (for  the  N unknowns  at  a 
point),  minimum-residual  Block  Gauss-Seidel,  and  Conjugate  Gradients.  The  user  specifies 
r,  the  number  of  multi-level  cycles,  and  m,  the  number  of  smoothing  iterations  in  each 
paxt  of  the  cycle. 

For  systems  of  PDEs,  the  linear  equation  matrices  are  nonsymmetric,  and  may  be  difficult 
to  solve;  there  is  no  guarantee  of  convergence.  Other  smoothing  methods  would  be  desir- 
able. Of  the  types  in  B2DE,  Gauss-Seidel  smoothing  has  been  most  satisfactory,  even  for 
problems  without  theoretical  reason  to  expect  convergence. 

Unless  u is  close  to  a solution  of  the  nonlinear  finite-element  equations,  u v may  be 
worse  thzin  u.  A damped  Newton’s  method  is  used:  replace  u by  u -t-  A v,  where  A is 
chosen  so  that  j]  r(u  + Av)  ||<||  r(u)  |1  where  ||  r ||  is  the  Euclideatn  norm  of  r.  This  avoids 
divergence  of  the  iterative  process,  but  does  not  guarantee  convergence. 

In  practice,  the  initial  triangulation  is  usually  too  coarse  to  give  the  desired  accuracy; 
local  refinement  is  needed  to  represent  the  solution  accurately  in  part  or  all  of  the  region. 
After  the  package  has  converged  to  an  approximate  solution  on  the  initial  triangulation, 
the  approximate  solution  may  be  used  to  estimate  the  error  and  produce  a new  and  finer 
triangulation  adaptively.  There  are  various  ways  to  estimate  the  error  in  each  triangle. 
For  the  present  work,  a local  method  is  iised. 

Consider  a single  triangle;  for  simplicity  we  suppose  its  vertices  are  1,  2,  and  3.  The 
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approximate  solution  inside  the  triangle  is 


3 

= X]  (11) 

m=l 

We  tentatively  divide  the  triangle  into  four  similar  triangles  by  connecting  the  midpoints 
of  the  three  sides.  On  the  four  triangles,  a potentially  more  accurate  u,-  is 

= XI  (12) 

m=l 

(The  primes  refer  to  the  finer  trizmgulation.)  We  keep  the  old  values  at  the  original  3 
vertices;  = aim,  ^ — 1»2,3.  Values  at  the  3 new  vertices  axe  foimd  by  solving  the 
linearized  finite  element  equations  for  m = 4,5,6.  The  error  in  the  triangle  can  be 
estimated  by  the  largest  change  at  any  of  the  points  4,  5,  and  6. 

This  estimate  is  not  accurate  in  absolute  value,  but  the  relative  estimates  in  different 
triangles  are  sufficiently  accurate  for  refining  the  triangulation. 

The  refined  triangulation  is  obtained  by  dividing  the  triangles  with  the  largest  estimated 
errors,  such  &s  all  triangles  with  estimated  error  larger  than  1/4  of  the  largest  error,  or  the 
worst  1/4  of  the  triangles.  If  the  refined  triangulation  adds  only  a few  more  triangles,  the 
user  may  choose  to  “compress,”  or  merge  the  new  triangulation  with  the  previous  one. 

For  the  refined  triangiilation,  the  previously  calculated  u is  probably  a good  guess  at  the 
new  u,  so  that  few  I'l’ewton  iterations  should  be  necessary.  The  hardest  work  is  usually 
done  on  the  initial  triangulation,  the  one  with  the  fewest  vertices. 

2.2  Program  Description 

B2DE  is  written  in  a private  version  of  Ratfor  [Kem75];  the  result  is  portable  Fortran  77. 
A few  routines  are  unavoidably  nonportable;  these  are  described  in  Section  5. 

The  PORT  library  [Fox78]  machine  constants  are  used  to  improve  portability.  All  work 
space  is  managed  using  the  PORT  library  stack. 

B2DE  may  be  viewed  as  a three-layer  program  package.  In  the  top  layer  are  the  programs 
written  by  the  user  to  describe  a particular  problem,  and  the  USRxx  programs  which  are 
user-modifiable.  In  the  middle  layer  Bite  programs  callable  by  users,  such  as  SOLxxx  and 
PLTxxx.  In  the  bottom  layer  are  programs  which  the  user  cannot  see.  The  user  does 
not  need  to  know  anything  about  the  bottom  layer  programs,  and  only  needs  to  know 
the  calling  interface  to  the  middle  layer  programs.  In  particular,  all  the  data  structure  of 
B2DE  is  invisible  to  the  user. 

B2DE  comes  with  interactive  driver  programs,  DRVxrx,  which  call  middle-  and  top-layer 
programs.  The  user  is  free  to  revise  the  driver  programs. 
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3.  Using  B2DE  with  Supplied  Drivers 

3.1  Specifying  Your  Problem 

To  iise  B2DE  on  your  program,  you  need  to  write  a main  program  and  subroutines  to 
evaluate  the  coefficients  o,  /,  and  g and  their  first  partial  derivatives.  The  easiest  way 
to  do  this  is  to  modify  2m  existing  program.  See  Appendix  A and  the  sample  programs 
on  the  distribution  tape  for  examples.  The  current  DRVxxx  driver  assumes  a data  file  to 
specify  the  geometry  of  the  region  and  a few  parameters.  See  Appendix  B and  the  sample 
problems  on  the  distribution  tape. 

Working  space  is  managed  inside  B2DE  using  the  PORT  stack  [Fox78].  Your  main  progrsm 
must  define  the  size  of  the  stack  and  initialize  it  using  iatkin,  as  in  Appendix  A.  The  size 
of  the  stack  determines  the  largest  problem  that  can  be  solved. 

Default  versions  of  the  USRxxx  routines  exist,  but  can  be  replaced  by  the  user.  See  Section 
3.4.  For  nonlinear  problems,  you  will  want  to  replace  USRGSS  to  provide  an  initiai  guess 
at  the  solution. 

Note  that  all  of  B2DE,  including  the  DRVxxx  routines,  uses  free-format  reading.  If  you 
make  a mistake,  such  as  putting  a decimal  point  in  an  integer,  or  putting  in  a random 
character,  you  will  be  informed  and  czm  make  a correction  if  you  are  running  interactively. 
If,  for  example,  a program  asks  for  four  integers  they  can  all  be  on  one  line,  or  on  separate 
lines.  Blank  lines  are  completely  ignored.  Input  lines  can  have  comments  in  them;  a # or 
a * begins  a comment,  which  continues  to  the  end  of  the  line. 

See  Section  6 for  the  examples. 

3.2  The  Interactive  Driver  - Finding  a Solution 

DRVxxx  programs  zu*e  drivers  that  make  calls  to  subroutines  in  B2DE  to  solve  the  PDFs. 
You  may  modify  or  replace  any  of  these  routines,  but  you  should  follow  the  model  of  the 
default  version  carefully. 

Your  main  program  should  be  modeled  on  USRMAI.  It  will  probably  call  DRVALL,  which  is 
the  driver  which  hemdles  the  input  phase  joid  then  calls  the  solution  driver.  The  default 
version  calls  USRINl,  then  SOLINP,  USRIN2,  USRIN3,  and  finally  DRVSOL.  DRVALL  assumes 
that  needed  files  have  been  opened  and  that  the  input/output  values  have  been  set  by  calls 
to  lOSET.  Note  that  ausr,  fusr,  and  gusr  must  be  declared  EXTERNAL;  the  names  can  be 
anything  you  like,  but  ausr,  fusr,  amd  gusr  will  be  used  in  this  manual. 

subroutine  DRVALL (ausr.  fusr.  gusr) 
external  ausr,  fusr.  gusr 

DRVSOL  is  the  driver  for  solving  the  PDEs.  DRVSOL  first  calls  TRMPLT  to  do  any  needed 
initializing  for  plotting,  then  calls  SOLREF  to  get  started.  The  main  part  of  DRVSOL  is  a 
loop  in  which  an  integer  is  read  and  one  of  nine  actions  is  taken.  Eight  of  the  actions 
are  calls  to  subroutines  DRVPRT,  DRVPLT,  SOLITR,  SOLREF,  SOLCMP,  SOLSAV,  SOLD  IF,  and 
SOLCHG.  All  th^e  are  described  separately,  and  may  be  called  by  the  iiser  in  a modified 
driver.  The  ninth  action  is  to  quit;  TRMPLT  is  then  called  to  do  any  needed  finishing  off  for 
plotting. 

subroutine  DRVSOL (ausr,  fusr,  gusr) 
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external  ausr,  fusr,  gusr 

The  remainder  of  this  section  is  a description  of  the  interactive  driver,  using  a run  from 
Example  1.  The  responses  from  B2DE  are  those  from  print  level  2. 

The  interactive  driver  is  menu-driven;  user  responses  are  one  or  more  numbers.  Except  at 
the  top  level,  entering  0 as  the  answer  to  a question  gets  you  to  the  next-higher  menu. 

On  entering,  B2DE  prints  the  level  (this  will  be  1 unless  an  old  solution  is  being  retrieved), 
the  time  taien  to  set  up  the  mesh  (grid) , the  number  of  vertices  in  the  level  1 mesh  (nv) , 
and  the  time  taken  to  do  initisd  set-up  (asmbnl). 

Then  the  main  menu  is  printed. 

refsol  entered  at  level  1 7 vertices 
grid  .033  sec.  16  vertices 
asmbnl  . 033  sec . 

l=print,  2*plot,  3=iter,  4=refine,  6=comp,  6=save,  7=diff,  8=end,  9=change 
The  possible  responses  will  be  described  separately. 

3.2.1  Printing 

The  DRVPRT  driver  program  is  called.  Its  main  printing  menu  DRVPRT  is 

l=8ummary,  2=timing,  3=storage,  4=vertices,  6=pde  variables,  6=settings 

Printing  is  to  the  terminal  unless  the  6=settings  menu  has  been  used  to  direct  printing 
and  plotting  to  files.  Each  of  the  responses  will  be  treated  in  order. 

3.2. 1.1  Summary 

PRTSUM  is  czdled.  For  each  level,  the  number  of  vertices  and  the  number  of  triangles  are 
printed.  For  each  PDE  and  each  level  for  which  refinement  has  been  attempted,  a measure 
of  the  solution  size,  unorm,  an  estimate  of  the  size  of  the  error,  enorm,  and  the  estimated 
number  of  correct  significant  digits  in  the  solution  are  printed.  The  error  and  the  estimated 
number  of  digits  aie  only  estimates;  they  are  not  guaranteed;  they  are  almost  certainly 
not  accurate. 

n Ivl  vertices  triangles  xmorm  enorm  digits 

11  62  114  7.04E-01  1.04E-02  1.83 

1 2 103  258 

3.2. 1.2.  Timing 

PRTTIM  is  called.  The  number  of  seconds  taken  in  each  of  the  major  parts  of  B2DE 
is  printed,  together  with  the  percentage  and  the  total,  assembly  is  the  time  taken  to 
calculate  Jacobian  matrices,  factor  is  the  time  taken  to  factor  the  level  1 matrices,  mg 
is  the  time  taken  to  solve  the  linear  equations  to  get  the  Newton  step,  rhs  check  is  the 
time  taken  to  evaluate  the  residuzd  to  see  if  a Newton  step  is  acceptable. 

timing  distribution 
seconds  percent 
plot  0.000  0.00 
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assembly 

2.200 

41.51 

factor 

0.210 

3.96 

mg 

0.180 

3.40 

rhs  check 

1.780 

33.58 

refine 

0.930 

17.55 

total 

5.300 

100.00 

3.2.1.3.  Storage 

PRTSTO  is  called.  The  maximum  number  of  triangles  and  vertices  are  printed,  together 
with  the  number  used  so  far.  The  maximum  number  depends  on  the  size  of  the  stack  used. 

103  vertices  used. 

2094  vertices  allowed. 

263  triangles  used. 

5583  tri2mgles  allowed. 

3.2. 1.4.  Vertices 


PRTVER  is  called.  The  coordinates  of  each  vertex  are  printed,  together  with  the  “parents” 
of  the  vertex;  the  parents  are  the  two  vertices  of  the  triangle  edge  which  was  bisected  to 
produce  the  vertex.  (The  original,  user-defined  vertices  are  their  own  parents.) 


parents 

X 

y 

1 

1 

1 

-5.00000E-01 

O.OOOOOE+00 

2 

2 

2 

O.OOOOOE+00 

5.00000E-01 

3 

3 

3 

5.00000E-01 

• 5.00000E-01 

4 

4 

4 

O.OOOOOE+00 

O.OOOOOE+00 

5 

5 

5 

5.00000E-01 

O.OOOOOE+00 

6 

6 

6 

O.OOOOOE+00 

-5.00000E-01 

7 

7 

7 

-5.00000E-01 

-5.00000E-01 

8 

2 

4 

O.OOOOOE+00 

2.50000E-01 

9 

1 

4 

-2.50000E-01 

O.OOOOOE+00 

10 

1 

2 

-3.53553E-01 

3.53553E-01 

11 

3 

4 

2.50000E-01 

2.50000E-01 

3.2. 1.5.  PDE  variables 

PDE  values  are  printed.  If  a transformation 

is  in  effect,  the  transformed  valueds  are  printed 

when  PDE  1 

is  asked  for. 

PRTPDE  is  called.  Its  main  menu  is  : 

l»on  grid,  2=at  vertices,  3=on  grid  as  scaled  integers 

Values  can  be  obtained  on  a regular  rectangular  grid  or  at  the  vertices.  Grid  values  can  be 
printed  as  floating-point  numbers  or  as  scaled  integers  (xiseful  for  quickly  getting  an  idea 
of  the  solution).  The  three  types  will  be  discussed  in  order. 

3.2.I.5.I.  On  Grid 

A PDE  is  evaluated  on  a rectangular  grid  aligned  with  the  i-  and  y-axes.  If  the  magni- 
fication is  1,  the  grid  is  on  a rectangle  defined  by  the  minimum  and  maximum  i and  y 
vertex  values.  If  the  magnifications  is  greater  than  1,  the  grid  covers  only  a portion  of  this 
rectangle. 
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The  next  menu  is: 

enter  pde  no.,  l=value  only,  2=gradient  only,  3-both 
Enter  the  PDE  number,  from  1 through  npde,  then  1,  2,  or  3. 
nx,  ny 

You  are  then  asked  for  the  number  of  grid  points  along  x ajid  along  y.  If  the  region 
is  not  rectangular,  or  if  its  edg^  are  not  aligned  with  the  axes,  some  grid  points  will 
fall  outside  the  region.  These  points  have  a large,  machine-dependent  number  printed  to 
denote  out-of-bounds;  in  these  examples,  the  number  is  1.7014lE-(-38. 

A typical  printing  is: 


X 

1 -5.00000E-01 

2 -2.50000E-01 

3 O.OOOOOE+00 

4 2.50000E-01 

5 5.00000E-01 

6 -6.00000E-01 


y 

-5.00000E-01 

-5.00000E-01 

-5.00000E-01 

-5.00000E-01 

-5.00000E-01 

-2.50000E-01 


u 

O.OOOOOE+00 

O.OOOOOE+00 

O.OOOOOE+00 

1.70141E+38 

1.70141E+38 

4.03474E-01 


3. 2. 1.5. 2.  At  Vertices 


You  axe  asked  which  PDE  to  evaluate.  Since  the  gradient  is  not  continuous  at  the  vertices, 
it  is  not  printed.  • 

enter  pde  no. 

A typical  printing  is: 


X 

1 -5.00000E-01 

2 O.OOOOOE+00 

3 5.00000E-01 

4 O.OOOOOE-^00 

5 6.00000E-01 

6 O.OOOOOE+00 


y 

O.OOOOOE+00 

5.00000E-01 

6.00000E-01 

O.OOOOOE+00 

O.OOOOOE+00 

-5.00000E-01 


u 

l.OOOOOE+00 

l.OOOOOE+00 

6.68726E-01 

7.04038E-01 

O.OOOOOE+00 

O.OOOOOE+00 


3. 2. 1.5. 3.  On  Grid  as  Sc^ed  Integers 

A PDE  is  evaluated  on  a rectangular  grid  aligned  with  the  x-  and  y-axes.  If  the  magni- 
fication is  1,  the  grid  is  on  a rectangle  defined  by  the  Tninirmim  and  maximum  x and  y 
vertex  vedues.  If  the  magnifications  is  greater  than  1,  the  grid  covers  only  a portion  of  this 
rectangle. 

The  next  menu  is: 

enter  pde  no.,  l=value  only,  2“gradient  only,  3=  both 
Enter  the  PDE  number,  from  1 through  npde,  then  1,  2,  or  3. 
nx,  ny 
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You  axe  then  asked  for  the  number  of  grid  points  along  x and  along  y.  If  the  region  is 
not  rectangular,  or  if  its  edges  axe  not  aligned  with  the  axes,  some  grid  points  will  fall 
outside  the  region.  These  points  have  9999  printed  to  denote  out-of-bounds.  Other  points 
are  scaled  to  lie  from  -100  to  100.  A typical  printing  is: 

X limits  -1.66667E-01  1.66667E-01 

y limits  -6.00000E-01  -1.66667E-01 

pde  1 

scaled  values  from  -6.085E-08  to  5.492E-01  scale  5.492E-01 


90 

90 

91 

93 

95 

97 

100 

78 

78 

80 

82 

85 

89 

92 

64 

65 

67 

70 

75 

80 

85 

49 

50 

54 

58 

65 

71 

77 

34 

35 

39 

44 

54 

62 

70 

17 

19 

21 

28 

42 

55 

65 

0 

0 

0 

0 

9999 

9999 

9999 

3.2. 1.6.  Settings 
The  next  menu  is 

l=color,  2=erase.  3=magnify,  4=window,  6*file  print /plot,  e^transform 

The  same  menu  is  available  from  the  plotting  menu.  All  six  menu  choices  afiFect  plotting, 
but  only  choices  3,  5,  and  6 affect  printing.  This  menu  is  discussed  in  detail  in  Section 
3.3.8,  under  plotting. 

3.2.2  Plotting 

Plotting  is  covered  in  great  detail  in  Section  3.3. 

3.2.3  Iterating 

Using  the  current  mesh,  do  one  or  more  Newton  iterations  on  the  nonlinear  finite  element 
equations  which  represent  the  PDEs.  In  each  iteration,  a Jacobism  matrix  emd  a right- 
hzmd-side  are  calculated;  the  solution  of  these  equations  is  the  Newton  correction.  If  the 
mesh  level  is  1,  the  linear  equations  are  solved  by  sparse  Gaussian  elimination.  If  the  mesh 
level  is  greater  than  1,  the  linear  equations  are  solved  by  a multi-level  iteration,  with  the 
level  1 equations  solved  by  spjirse  Gaussian  elimination. 

If  the  PDEs  are  linear,  one  iteration  is  suflBcient  if  the  linear  equations  are  solved  sufficiently 
accurately.  The  prompt  message  is 

enter  niter,  ieps,  relerr.  abserr.  jsmax.  param 

niter  is  an  integer,  the  maximum  number  of  Newton  iterations  to  be  done.  Fewer 
will  be  done  if  the  accuracy  requirements  are  met. 

ieps  is  an  integer  which  governs  the  accuracy  with  which  the  linear  equations  arc 
solved.  Let  ||  r H be  the  Euclidean  norm  of  the  initial  right-hand-side;  this  is  also 
the  residual  of  the  nonlinear  equations  and  the  initial  residual  of  the  linear  equations. 
The  multi-level  iterations  are  stopped  when  the  new  residual  is  ||  r j|  or  less,  or 

when  maxr  multi-level  iterations  have  been  done,  maxr  is  set  by  subroutine  USRIN3. 
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Tlie  value  of  ieps  is  immaterial  if  the  level  is  1,  since  level  1 equations  axe  solved  by 
sparse  Gaussian  elimination. 

If  ieps  is  too  small,  the  linear  equations  will  not  be  solved  accurately  enough,  and  more 
iterations  may  be  needed.  If  ieps  is  too  large,  more  work  may  be  used  in  solving  the 
linear  equations  than  is  justified  by  the  accuracy  of  the  approximate  solution,  relerr 
and  abserr  are  real  numbers  which  govern  the  accuracy  with  which  the  nonlinear 
equations  are  solved.  Newton  iterations  are  stopped  when  the  residual  is  small  enough, 
when  II  r 1|<=  (relerr  ||  u ||  -fabserr),  where  ||  u ||  is  the  Euclidean  norm  of  the 
current  approximate  solution. 

j smax  is  an  integer  which  governs  the  damped  Newton  solution.  If  v is  the  Newton 
correction,  the  new  approximate  solution  is  u + tv.  First  t = to  is  tried,  where  to  is 
the  value  returned  by  USRSTP,  then  to/2,  then  to/4,  and  so  on,  until  the  new  residual 
is  small  enough;  that  is,  imtil  ||  r(u-f-tv)  |1<||  r(u)  ||.  A total  of  jsmax  tries  are  made. 

param  is  a real  parameter  for  the  convenience  of  the  user;  it  is  put  into  a common 
block  and  can  be  accessed  by  the  user.  It  may  be  useful  in,  for  example,  USRREG.  The 
common  block  is 

COMMON  /USRCOM/  PARAM 

If  the  nonlinear  iterations  have  not  converged  by  the  niter-th  iteration,  DRVSEE  is  called 
to  give  the  user  an  opportimity  to  look  at  the  Newton  correction  before  the  correction  has 
been  trated.  The  menu  is 

l=print,  2=plot,  3=continue 

Menu  choices  1 and  2 are  exactly  the  same  as  choices  1 amd  2 of  the  main  menu,  except  that 
the  user  prints  or  plots  the  current  Newton  correction  rather  than  the  current  approximate 
solution.  The  current  magnification,  window,  and  color  are  used.  Use  menu  choice  3 to  go 
on  to  test  the  Newton  correction. 

B2DE  has  provision  for  user-supplied  adjustments.  For  example,  the  PDEs  may  depend 
on  the  solution  to  some  nonlinear  equation  formally  separate  from  the  PDEs.  After  u has 
been  replaced  by  u -t-  tv,  but  before  the  residual  is  calculated,  USRADl  is  called.  After  a 
satisfactory  value  of  t has  been  found,  and  after  the  residual  is  calculated,  USRAD2  is  called. 

3.2.4  Refining 

Using  the  current  approximate  solution,  estimate  the  error  in  the  solution  and  adaptively 
refine  the  mesh.  B2DE  first  displays  the  current  level  number,  the  number  of  vertices,  and 
a crude  guess  at  the  maximum  error  in  the  solution  for  each  PDE.  A typical  example  is: 

refsol  re-entered  at  level  2 26  vertices 

asmbnl  . 050  sec . 

1 estimated  maximum  error  8.01E-03 

B2DE  then  gives  you  an  opportimity  to  plot  the  estimated  error  for  each  PDE.  When 
you  are  finished,  if  npde  is  2 or  more,  B2DE  gives  you  an  opportunity  to  plot  the  total 
estimated  error  for  all  the  PDEs  combined. 

B2DE  then  shows  the  estimated  Euclidean  error  for  each  PDE,  shows  how  many  triangles 
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will  be  refined  because  the  local  estimated  error  is  high  enough,  and  then  refines.  The 
estimated  error  and  the  number  of  correct  digits  axe  almost  certainly  wrong,  but  are 
indicative  of  the  progress  of  the  solution  from  level  to  level. 

If  any  of  the  weights  returned  by  USRWTS  is  negative,  B2DE  asks  for  the  weights  to  use 
and  asks  if  you  want  to  refine. 

level  1 pde  1 unorm,  enorm  7.06E-01  2.86E-02  eat.  digits  1.39 

min  error  / mzoc  error  1.44  percent 

4 of  30  triangles  ( 13.33)  above  60.00  percent 

9 of  30  triangles  ( 30.00)  above  25.00  percent 

8 of  30  triangles  ( 26.67)  above  37.50  percent 

B2DE  tells  you  the  number  of  vertices  at  the  new  level. 

refine  .700  sec.  level  2 52  vertices 

3.2.5  Compressing 

The  vertices  of  the  highest-level  mesh  are  merged  with  those  of  the  next-to-highest- level 
mesh,  and  the  level  number  is  reduced  by  1. 

This  option  is  useful  in  cases  where  a refinement  generates  only  a few  new  vertices,  as  may 
happen  in  problems  with  a singularity  in  the  solution. 

3.2.6  Saving  the  solution 

The  current  approximate  solution  is  saved  on  a file.  The  file  must  have  been  opened  by 
the  user  before  calling  B2DE.  Let  n = ioget(8). 

If  n > 0,  the  information  is  written  as  unformatted  (binary)  information  on  logical  file  \mit 
n.  If  n<  0,  the  information  is  saved  as  formatted  information  on  logical  file  unit  -n.  The 
prompt  message  is 

enter  the  number  of  significant  figures  to  save 

After  the  user  enters  the  number  of  significant  figures  to  save,  B2DE  constructs  an  appro- 
priate format  and  writes  the  file. 

Formatted  files  tadte  more  space  to  write,  but  can  be  used  to  transport  data  from  one 
computer  to  another. 

3.2.7  Looking  at  the  Difference 

Look  at  the  change  in  the  current  solution  since  the  last  time  the  mesh  was  refined,  or 
since  the  start  of  the  current  session  with  B2DE,  whichever  is  later.  DRVSEE  is  called. 

The  menu  is 

l=print , 2=plot , 3*continue 

Menu  choices  1 and  2 are  exactly  the  same  as  choices  1 aind  2 of  the  main  menu,  except 
that  the  user  can  print  or  plot  the  ciirrent  difference  rather  than  the  ctirrent  approximate 
solution.  Use  menu  choice  3 to  go  back  to  the  main  menu. 

3.2.8  Quitting 
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Quit. 

3.2.9  Changing  Parameters 

Change  user-defined  parameters,  internal  B2DE  parameters,  or  change  which  PDE  vari- 
ables are  being  iterated  on. 

First  USRCHG  is  called;  in  USRCHG,  the  user  may  change  any  xiser-defined  parameters  as 
desired. 

If  npde  is  2 or  more,  the  user  may  keep  the  values  of  one  or  more  of  the  PDE  variables 
fixed,  and  calculate  Newton  corrections  for  the  remaining  PDE  variables.  If  npde  is  2 or 
more,  B2DE  next  prints 

PDE  variables  - 0 for  fixed,  1 for  varying 
2nd  prints  the  current  status  of  each  of  the  npde  PDE  variables, 
enter  new  values 

The  user  must  enter  the  new  values  as  integers:  1 means  the  PDE  is  to  vary,  and  einy  other 
integer  means  the  PDE  is  held  fixed.  At  least  one  1 must  be  entered. 

The  next  prompt  is 

enter  1 to  change  pareimeters , 0 to  return 
If  1 is  entered,  the  next  prompt  is 

parameters:  maxm,  maxr,  ittype,  ithrsh,  iprsw 

and  the  current  values  are  printed  out.  Enter  the  new  values,  all  integers.  Refer  to  Section 
3.4.9,  USRIN3,  for  the  meanings  of  the  parameters.  m2Lxm,  maxr,  and  ittype  all  govern 
multi-level  iteration,  ithrsh  governs  refinement,  iprsw  governs  the  amoimt  of  printing  as 
the  solution  progresses. 

3.3  The  Interactive  Driver  - Graphics 

Seeing  pictures  is  necessary  if  the  user  is  (l)  to  understand  what  is  happening  during  the 
iterative  process  leading  to  a solution  and  (2)  to  understand  the  properties  of  the  solution. 

B2DE  gives  the  user  ample  opportunity  to  msike  plots.  There  are  standard  plots  available 
of  several  kinds.  The  variables  plotted  axe  either  a PDE  variable  or  the  magnitude  of  its 
gradient.  Since  the  author  cannot  guess  everything  the  user  might  want  to  plot,  B2DE 
allows  the  user  to  define  problem-dependent  plot  variables  (Section  3.4.12).  Plots  may 
be  made  to  the  user’s  graphics  terminal,  or  the  vectors  may  be  written  to  a file  for  later 
processing. 

Most  of  the  graphics  progrsim  code  is  independent  of  the  computer  used  and  of  the  graphics 
output  device.  The  dependent  parts  are  treated  in  the  installation  part  of  the  manual. 
Section  5. 

Five  kinds  of  plots  are  available:  triangles,  contours,  surfaces,  profiles,  and  flow  lines.  Each 
plot  type  is  done  by  a subroutine  which  may  be  called  directly  from  a xiser’s  program.  Each 
plot  type  also  has  an  interactive  driver  which  may  be  called  directly  from  a user’s  program. 
There  is  &n  overall  interactive  plot  driver. 
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The  interactive  driver  is  DRVPLT.  The  user  may  wish  to  modify  or  replace  this  progrsim. 
DRVPLT  asks  interactively  for  user  input,  then  calls  any  or  all  of  DRVTRI,  DRVCON,  DRV- 
SUR,  DRVPRO,  DRVFLW,  and  DRVSET.  The  user  may  wish  to  modify  or  replace  2iny  of  these 
programs. 

Most  of  the  DRVxxx  programs  call  the  corresponding  PLTxxx  programs:  PLTTRI,  PLTCON, 
PLTSUR,  PLTPRO,  and  PLTFLW.  The  calling  sequences  for  the  PLTxxx  programs  are  described 
in  Section  4.3;  PLTxxx  programs  may  be  called  directly  by  the  user,  but  are  not  meant 
to  be  changed  by  the  user.  Some  PLTxxx  programs  require  work  space,  which  is  obtained 
from  the  stack;  if  insufficient  stack  space  is  available,  the  plot  is  skipped. 

DRVSET  calls  any  of  six  sub-driver  programs:  DRVCOL,  DRVERA,  DRYMAG,  DRVWDW,  DRVFPL,  and 
DRVSHF.  The  DRVxxx  programs  call  the  corresponding  SETxxx  programs:  SETERA,  SETCOL, 
SETMAG,  SETWDW,  SETFPL,  and  SETSHF.  The  calling  sequences  for  the  SETxxx  programs  are 
described  in  Section  4.1;  SETxxx  progrzuns  may  be  called  directly  by  the  user,  but  are  not 
meant  to  be  changed  by  the  user. 

Plots  can  be  drawn  to  the  full  screen  or  plotting  surface,  or  can  be  drawn  in  a subwindow. 
Contour  plots,  flow  line  plots,  and  triangle  plots  can  be  magnified. 

All  plots  axe  of  vectors  only,  with  no  plot  labeling.  If  the  hardware  \ised  supports  it,  erasing 
between  plots  is  optional.  If  the  hardware  used  supports  it,  the  line  type  or  color  of  the 
plots  can  be  chamged. 

The  following  description  of  plotting  assumes  the  overall  interactive  driver,  DRVPLT.  All 
the  DRVxxx  programs  are  menu-driven  and  us*e  a free-format  input  package.  As  in  the  rest 
of  B2DE’s  interactive  drivers,  responses  of  0 return  the  user  to  the  next  higher-up  menu. 
Erroneous  responses  are  dealt  with  reasonably;  error  messages  should  be  self-explanatory. 
The  top-level  menu  of  DRVPLT  is: 

l^triangles,  2=contour3,  3=surface,  4=profilea,  6=flow  lines,  6=3ettings 

The  response  should  be  a non-negative  integer.  A response  of  0 meems  to  go  to  the  next 
higher  menu.  Other  responses  outside  the  indicated  range  are  ignored.  The  responses  will 
be  treated  in  order. 

3.3.1  Triangle  plots 

The  DRVTRI  driver  is  called.  The  next  menu  is  : 
enter  level 

Any  existing  level  of  the  triangulation  can  be  plotted.  The  current  magnification  and 
window  are  used.  The  plot  is  in  the  current  color. 

After  the  plot  is  done,  the  program  may  wait  for  the  user’s  signal.  The  signal,  if  any,  is 
that  required  by  the  terminal-dependent  program,  PLTUTL.  After  the  signal,  the  time  taken 
is  displayed,  and  the 

enter  level 

menu  comes  up  again.  To  get  back  to  the  previous  menu,  enter 
0 0 
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3.3.2  Contour  Plots 

The  DRVCON  driver  is  called.  The  next  menu  is  : 

pde  no.,  l=f;mction/2=ab8(grad) , l=fill/2=no  fill 

Enter  the  PDE  number,  from  1 through  iV,  the  number  of  PDEs  in  the  problem.  Enter  1 
for  contours  of  the  PDE  value,  2 for  the  absolute  value  of  its  gradient.  Enter  1 to  fill  in 
the  contours  with  different  colors  in  the  current  palette,  or  2 to  just  draw  the  contours  in 
the  current  line  color.  The  current  magnification  and  window  are  used.  If  fill  is  chosen, 
the  next  menu  is  : 

start  with  color  number 

Enter  the  first  color  of  the  color  palette,  which  may  be  from  0 on  up.  The  space  between 
the  bottom  two  contours  will  be  filled  with  this  color;  succeeding  spaces  will  be  filled  with 
successive  colors.  The  next  menu  is  : 

number  of  contours 

The  number  of  contours  is  limited  to  101  by  PARAMETER  (MAXCON  = 101)  in  program 
DRVCON.  The  magnification  is  treated  as  for  triangle  plots. 

If  nc,  the  number  of  contours,  is  positive,  the  program  finds  the  maximum  and  minimum 
values  of  the  plot  variable  and  spaces  the  contours  equally  between,  but  not  including,  the 
extreme  values.  Both  the  extreme  values  and  the  contour  values  are  displayed. 

If  nc  is  0,  the  extreme  values  are  displayed,  but  no  plot  is  done. 

If  nc  is  negative,  the  user  is  prompted  for  extreme  values  for  contours.  There  are  -nc 
contours,  equally  spaced  between  and  including  the  extreme  values. 

If  no  fill  is  chosen,  the  next  prompt  is  : 

number  of  contours 

and  the  replies  are  as  above. 

After  the  plot  is  done,  the  time  taken  is  displayed,  and  the 
pde  no.,  l»function/2»abs(grad) , l=fill/2=no  fill 
menu  comes  up  again.  To  get  back  to  the  previoiis  menu,  enter 
0 0 0 

3.3.3  Surface  Plots 

The  DRVSUR  driver  is  called.  The  current  color  and  window  are  used.  The  magnification 
is  1.  The  next  menu  is  : 

pde  no.,  l=function/2=abs(grad) 

Enter  the  PDE  number,  from  1 through  N.  Enter  1 for  a plot  of  the  PDE  value,  2 for  a 
plot  of  the  absolute  value  of  its  gradient.  The  next  menu  is  : 

level,  l=surface  triemgles,  2=surface  contours,  3=surface  grids 

Enter  the  triangulation  level  desired  for  the  plot,  and  the  type  of  stirface  plot  desired.  The 
absolute  value  of  the  level  is  used.  In  all  three  types,  the  function  is  plotted  as  a parallel 
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projection  plot  (no  perspective);  the  height  of  the  surface  above  or  below  the  x-y  plane 
represents  the  function  value.  The  surface  can  be  viewed  from  any  angle.  Hidden  lines 
are  omitted  unless  the  level  is  negative;  then  hidden  lines  are  plotted.  Plotting  without 
removing  hidden  lines  requires  much  less  computation,  and  can  be  much  faster. 

The  three  types  will  be  discussed  in  order. 

3.3. 3.1  Surface  triangles 

The  surface  is  defined  by  the  triangles  of  the  specified  level.  The  next  menu  is  : 
enter  viewing  direction  — x,  y,  z 

The  plot  is  displayed  as  if  seen  by  a viewer  looking  towards  the  surface  along  a line 
extending  from  the  origin  through  the  point  (x,  y,  2).  After  the  plot  is  done,  the  time 
taken  is  displayed,  and  the 

pde  no.,  l=f xmct ion/ 2=abs (grad) 

menu  comes  up  again.  To  get  back  to  the  previous  menu,  enter 

0 0 

3.3.3.2  Surface  contours 

The  surface  is  defined  by  the  contour  lines  drawn  on  triangles  of  the  specified  level.  Triangle 
boundaries  are  not  shown.  The  first  menu  is  again  : 

pde  no.,  l=function/2=ab8(grad) 

and  the  responses  are  as  before.  The  next  menu  is  : 

number  of  conto\irs 

Responses  to  this  menu  are  like  those  for  the  planar  contour  plot,  except  no  magnification  is 
allowed.  The  number  of  contours  is  limited  to  202  by  PARAMETER  (MAX  » 202)  in  program 
DRVSUR. 

The  next  menu  is  : 

enter  viewing  direction  — x,  y,  z 

Responses  are  as  before.  After  the  plot  is  done,  the  time  tadcen  is  displayed,  and  the 
pde  no.,  l=function/2=abs(grad) 
menu  comes  up  again.  To  get  back  to  the  previous  menu,  enter 
0 0 

3.3.3.3  Surface  grids 

The  surface  is  defined  by  grid  lines  parallel  to  the  x zind  y aoces,  drawn  on  triangles  of  the 
specified  level.  Triangle  boimdaries  Me  not  shown.  The  first  menu  is  again 

pde  no.,  l=function/2=abs(grad) 

and  the  responses  are  as  before.  The  next  menu  is 

number  of  x lines,  n\imber  of  y lines 
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The  X lines  £ire  equally  spax:ed  along  the  x-axis,  from  the  minimum  x- value  to  the  maximum 
X- value;  the  y lines  are  similarly  spaced  along  the  y-axis.  The  total  number  of  grid  lines,  x 
plus  y,  is  limited  to  202  by  PARAMETER  (MAX  = 202)  in  program  DRVSUR.  The  next  menu 
is: 

enter  viewing  direction  — x,  y,  z 

Responses  are  as  before.  After  the  plot  is  done,  the  time  taken  is  displayed,  and  the 
pde  no.,  l=function/2=abs(grad) 
menu  comes  up  again.  To  get  ba^k  to  the  previous  menu,  enter 
0 0 

3.3.4  Profile  Plots 

The  DRVPRO  driver  is  called.  The  current  color  zmd  window  are  \ised.  The  next  menu  is  : 
enter  pde  no. , number  of  profiles 

Enter  the  PDE  number,  from  1 through  iV,  and  the  number  of  profiles  desired.  A profile 
is  a slice  through  a surface  plot,  viewed  from  the  side.  The  slices  Eire  parallel  and  equally 
spaced.  The  next  menu  is  : 

enter  indices  of  profile  plane (s)  x,  y 

The  slices  are  parallel  to  a line  running  from  the  origin  to  the  point  (x,  y).  Viewing  is 
from  the  direction  as  defined  for  surface  plots,  (y,  -x,  0),  perpendicular  to  the  slices.  After 
the  plot  is  done,  the  time  taken  is  displayed,  and  the 

enter  pde  no.,  number  of  profiles 

menu  comes  up  again.  To  get  back  to  the  previous  menu,  enter 

0 0 

3.3.5  Flow  Line  Plots 

The  DRVFLW  driver  is  called.  Flow  line  plots  attempt  to  describe  the  fiux,  or  vector  field, 
a,  Vui.  The  user  specifies  a line  segment,  or  set  of  connected  segments;  the  positive  and 
negative  fiuxes  across  the  line  are  calculated;  points  along  the  line  are  spaced  evenly  in 
the  fiux;  and  flow  lines  2U“e  drawn  both  ways  from  these  points.  Flow  lines  axe  terminated 
at  boimdaries  of  the  region  or  when  the  fiow  lines  would  re-enter  a triangle.  The  current 
magnification  and  window  are  used.  The  plot  is  in  the  current  color. 

Since  the  solutions  are  composed  of  linear  finite  elements  on  triangles,  Vu»  is  constant 
within  each  triangle.  Gradients  tend  not  to  go  smoothly  to  zero  when  they  should,  but 
overshoot  instead;  this  makes  the  fiow  lines  stop. 

The  next  menu  is  : 

enter  pde  no.,  number  of  points,  number  of  flow  lines 

The  PDE  number  can  be  from  1 to  N.  The  number  of  points  is  the  number  in  the  set  of 
connected  line  segments,  two  or  more.  The  next  menu  is  : 

enter  # x-values  followed  by  # y-values 
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where  # is  the  number  of  flow  lines  just  entered.  The  total  flow  in  a positive  and  in  a 
negative  directions  across  the  broken  line  is  printed,  and  then  the  flow  lines  are  drawn. 
After  the  plot  b done,  the  time  taken  is  displayed,  and  the 

enter  pde  no.,  number  of  points,  number  of  flow  lines 

menu  comes  up  again.  To  get  back  to  the  previous  menu,  enter 

0 0 0 

3.3.6  Settings 

The  DRVSET  driver  is  called.  Several  values  can  be  set;  these  apply  to  future  plots  until 
the  settings  are  changed.  The  menu  is 

l=color.  2=erase.  3=magnify.  4=window.  5=file  print/plot,  6=transform 
Setting  the  Color 

The  DRVCOL  driver  is  called.  Some  terminals  can  support  vectors  drawn  in  different  colors 
or  line  types.  The  menu  is  : 

l^change  color  number.  2=:change  color  palette 

If  1 is  entered,  the  next  menu  is: 

enter  new  color  number 

SETCOL  is  called,  which  calls  TRMCOL,  and  then  the  settings  menu  appears. 

If  2 is  entered  in  response  to 

l*chaoige  color  number,  2=change  color  palette 
the  system-dependent  subroutine  TRMPAL  is  called.  Then  the  settings  menu  appears. 
Setting  the  Erase  Flag  and  Erasing  the  Screen 

The  DRVERA  driver  is  called.  Some  terminals  support  optionad  erasing  between  plots.  If 
the  erase  flag  is  ON,  the  screen  is  erased  between  plots;  if  the  erase  flag  is  OFF,  the  screen 
is  not  erased  between  plots.  You  might  want  to  plot,  for  instance,  s\irface  contours  on  the 
same  plot  as  a surface  grid,  or  planar  contours  of  more  than  one  PDE.  The  menu  is: 

l=set  erase  flag  off.  2=set  erase  flag  on.  3=erase  screen 

If  option  3 is  chosen,  the  screen  is  erased  immediately,  but  the  erase  flag  is  not  changed. 
There  are  no  lower  menus.  After  the  desired  option  is  chosen,  SETERA  is  called  and  the 
settings  menu  appears. 

Setting  the  Magnification 

The  DRVMAG  driver  is  called.  B2DE  can  make  magnified  triangle,  flow  line,  and  contour 
plots.  The  parameters  necesary  are  the  magnification  and  the  point  which  will  be  the 
center  of  the  magnified  plot.  The  magnification  also  affects  printing  of  PDE  variables  on 
a grid,  either  as  real  values  or  as  scaled  integers.  The  initial  menu  is: 

enter  magnification 

If  the  magnification  entered  is  1 or  less,  the  magnification  is  set  to  1 and  you  are  returned 
to  the  main  plot  driver  menu.  If  the  magnification  is  greater  than  1,  the  next  menu  is: 
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enter  x-center  am.d  y-center  (0  to  100) 

An  unmagnified  plot  fits  in  a rectangle  aligned  with  the  x-  and  y-axes,  defined  by  the 
TniniTTmm  £ind  maximum  x and  y vertex  values.  For  magnifications  greater  than  1,  the 
center  of  the  magnified  plot  is  defined  with  respect  to  this  rectangle;  the  center  is  defined  by 
the  percentage  distance  from  the  minimum  to  the  maximum  coordinates  of  the  rectangle. 
For  plotting,  the  center  may  be  adjusted  to  fill  as  much  of  the  screen  as  possible.  For 
example,  if  the  center  is  given  as  100  100,  the  upper  right  comer  of  the  rectangle  will  be 
plottd  to  the  upper  right  comer  of  the  screen  rather  than  to  the  center  of  the  screen. 

SETMAG  is  called  and  the  settings  menu  appears. 

Setting  the  Window 

The  DRVWDW  driver  is  called.  B2DE  can  draw  its  plots  in  the  entire  plotting  surface  as 
defined  by  TRMSCL,  or  can  use  a subwindow.  The  prompt  is 

enter  x-left,  x-right,  y-bottom,  y-top  (from  0 to  100) 

The  numbers  entered  are  interpreted  aa  percentages  of  the  full  window.  For  exaonple, 
0 50  50  100  means  the  upper-left  quarter  of  the  full  window.  SETWDW  is  called  and  the 
settings  menu  appears. 

Printing  and  Plotting  to  Files 

The  DRVFPL  driver  is  called.  Printing  and  plotting  can  be  done  to  the  screen  of  the  terminal 
or  to  a file.  One  fiag  governs  both  printing  and  plotting.  The  prompt  is 

l=plot/print  to  screen,  2=plot/print  to  file 

SETFPL  is  called  and  the  settings  menu  appears. 

Users  may  write  their  own  programs  to  read  the  plot  file  and  produce  plots  on  other 
devices.  The  format  of  the  file  is  as  follows.  The  first  line  is  produced  by  the  terminal- 
dependent  program  PLTUTL,  and  may  be  changed  by  the  user.  The  default  version  writes 
a line  containing  the  minimum  and  maximum  integer  values  produced  by  the  plot, 

xmin  ymin  xmax  ymax 

in  (417)  format.  These  values  come  from  the  terminal-dependent  program  TRMSCL,  eind 
may  be  changed  by  the  user;  the  default  values  are  500,  9500,  500,  and  7250. 

The  plot  vectors  are  next  in  the  file.  Each  is  on  a separate  line, 

xl  yl  x2  y2 

in  (417)  format.  The  end  of  the  plot  is  marked  by  a plot  vector 
-1  -1  “1  -1 

which  is  produced  by  TRMPLT,  and  which  may  be  changed  by  the  user. 

Transforming  a Solution 

The  DRVSHF  driver  is  called.  The  user  may  specify  transformations  of  PDEs  and  print  or 
plot  the  transformed  values.  The  prompt  is 

enter  index  of  transformation 
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SETSHF  is  called,  which  calls  USRSHF,  and  the  settings  menu  appears. 

3.4  Easy  Modifications  - USlbooc  Programs 

Default  routine  USRxxx  aire  supplied;  these  routines  are  designed  to  be  modifiable  by  the 
user.  They  provide  ezisy  ways  to  alter  the  performance  of  B2DE  without  going  inside 
B2DE.  Users  should  read  the  default  versions  carefully  before  modifying  the  programs. 

3.4.1  USRMAI  is  the  Fortran  “main”  program.  USRMAI  declares  storage,  sets  up  for  input 
and  output  (I/O),  and  calls  a driver  routine.  The  a,  /,  and  g programs  should  be  declared 
EXTERNAL  in  this  program.  The  stack  must  be  declared  and  initialized  in  the  main  program. 
1/ O values  must  be  set  by  calling  lOSET  for  all  nine  values.  These  values  may  be  retrieved 
by  calling  lOGET;  B2DE  routines  needing  to  do  I/O  get  file  numbers  in  this  way. 

io(l)  is  the  file  number  to  be  used  for  reading  from  the  terminal. 

io(2)  is  the  file  number  to  be  used  for  writing  to  the  terminal. 

io(3)  is  the  terminal  type;  it  is  mainly  used  by  graphics  programs.  If  io(3)  is  non- 
positive, the  session  is  assumed  to  be  batch,  not  interactive.  All  plots  will  be  made 
to  files,  and  FREEIN,  the  free-format  input  routine,  will  abort  if  am  input  line  is  bad. 

io(4)  is  the  file  number  to  be  used  for  reading  input  data;  it  is  only  used  in  USRxxx 
programs. 

io(5)  is  the  file  number  to  be  used  for  writing  plot  data,  if  any  platting  is  done  to 
files  rather  than  to  the  terminal. 

io(8)  is  the  file  number  to  be  used  for  writing  print  data,  if  any  printing  is  done  to 
files  rather  than  to  the  terminal. 

io(7)  is  the  file  number  to  be  used  for  reading  the  old  solution  data,  if  an  old  solution 
is  to  be  read.  If  io(7)  is  positive,  the  solution  data  is  assumed  to  have  been  written 
as  unformatted  data.  If  io(7)  is  negative,  the  solution  data  is  assumed  to  have  been 
written  as  formatted  data. 

io(8)  is  the  file  number  to  be  used  for  writing  the  new  solution  data,  if  a solution 
is  to  be  saved.  If  io(8)  is  positive,  the  solution  data  will  be  written  as  unformatted 
data.  If  io(8)  is  negative,  the  solution  data  will  be  written  as  formatted  data. 

io(9)  is  the  file  number  to  be  used  by  B2DE  for  writing  scratch  data.  Data  is  written 
to  and  read  from  io(9)  as  unformatted  data. 

Data  files  corresponding  to  file  numbers  io(4),  io(5),  io(6),  io(7),  and  io(8)  should 
be  opened  if  there  is  any  chance  they  will  be  used.  File  io(9)  is  always  used  and  should 
be  opened. 

DRV  ALL  or  another  driver  should  be  called. 

Any  files  opened  should  be  closed. 

3.4.2  USRADl  is  called  to  make  ainy  needed  adjustments  after  a tentative  Newton  step  has 
been  taken,  but  before  the  new  residual  is  calculated.  The  default  version  does  nothing. 

subroutine  USRADl 
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3.4.3  USRAD2  is  called  to  make  any  needed  adjustments  after  a Newton  step  has  been 
accepted,  after  the  new  residual  is  calculated.  The  default  version  does  nothing. 

subroutine  USRAD2 

3.4.4  USRCHG  is  called  when  the  9=change  selection  from  the  main  menu  is  taken.  It 
should  return  .TRUE,  if  the  changes  made  affect  the  residuals  of  the  PDEs  so  they  must 
be  recalculated;  it  should  return  .FALSE,  if  the  residuals  need  not  be  recalculated.  The 
default  version  returns  . FALSE . 

logical  function  USRCHG (junk) 

junk  is  a dummy  parameter. 

3.4.5  USRGET  is  called  when  an  old  solution  is  being  read  in,  to  read  anything  that  was 
saved  by  USRSAV.  The  default  version  does  nothing.  USRGET  should  read  values  in  the 
same  order  and  with  the  same  formats  as  USRSAV  writes.  USRGET  is  always  called  with 
isin  equal  to  io(7). 

subroutine  USRGET (isin) 
integer  isin 

isin  is  an  input,  the  hie  number  to  read  from. 

3.4.6  USRGSS  is  called  to  get  the  initial  guess  at  the  solution.  The  default  version  guesses 
0;  this  is  adequate  for  linear  problems,  but  is  probably  not  adequate  for  most  nonlinear 
problems.  USRGSS  Is  C2klled  for  each  point, in  the  initial  mesh. 

subroutine  USRGSS (npde,  x,  y,  i.  u) 
integer  npde.  i 
real  x,  y,  u(npde) 

npde  is  the  number  of  PDEs. 

X and  y are  the  coordinates  of  the  point. 

i is  the  number  of  the  user-supplied  triangle  in  which  the  point  lies, 
u is  an  output,  the  guess  at  the  solution  at  the  point. 

3.4.7  USRINl  is  the  first  input  routine.  It  returns  six  values.  The  default  version  calls 
I0G£T(4)  to  get  a file  number,  and  calls  FREEIN  to  read  the  values  from  the  file.  USRINl 
is  called  by  DRV  ALL;  if  DRVALL  is  replaced  by  the  user,  USRINl  may  not  be  needed. 

subroutine  USRINl (oldsol,  nt.  nv.  nc,  nb.  npde) 
logical  oldsol 

integer  nt , nv . nc , nb . npde 

oldsol  should  be  .TRUE,  if  an  old  solution  is  to  be  read  in,  and  .FALSE,  if  a new 
solution  is  to  be  started  from  the  beginning.  If  oldsol  is  . TRUE . , the  remaining  five 
values  do  not  need  to  be  returned;  the  old  values  sure  used  by  B2DE. 

nt  is  the  number  of  triangles  in  the  user-supplied  coarsest  mesh. 

nv  is  the  number  of  vertices  in  the  user-supplied  coarsest  mesh. 

nc  is  the  number  of  curved  boundary  edges  in  the  user-supplied  coarsest  mesh. 
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nb  is  the  number  of  boundary  edges  in  the  user-supplied  coarsest  mesh, 
npde  is  the  number  of  PDEs. 

3.4.8  USRIN2  is  the  second  input  routine.  It  is  called  to  obtain  the  geometry  of  the  user- 
supplied  coarsest  mesh.  It  is  called  only  if  oldsol  is  . FALSE . . The  default  version  calls 
lOGET  (4)  to  get  a file  number,  and  calls  FREEIN  to  read  the  values  from  the  file. 

subroutine  USRIN2(nt,  nv,  nc,  nb,  npde,  Iref,  vxin,  vyin, 

* xmid,  ymid,  itnode,  ibdry) 
integer  nt,  nv,  nc,  nb,  npde 

real  vxin(nv) , vyin(nv) , xmid(nc) , ymid(nc) 
integer  ibdryCnb,  npde+4) , itnodeO,  nt) , Iref(nv) 

nt,  nv,  nc,  nb,  and  npde  are  all  inputs;  they  are  as  described  under  USRINl.  The 
remainder  of  the  arguments  are  outputs. 

Iref  (j)  is  the  initial  refinement  level  for  vertex  j. 

vxin(j)  is  the  x-coordinate  of  vertex  j. 

vyin(j)  is  the  y-coordinate  for  vertex  j. 

xmid(j)  is  the  x-coordinate  for  the  midpoint  of  curved  edge  j,  if  nc  > 0. 
ymid(j)  is  the  y-coordinate  for  the  midpoint  of  curved  edge  j , if  nc  >0. 
itnode (l,k),  itnode (2, k),  and  itnode (3, k)  are  the  three  vertices  for  triangle  k. 
ibdry  (m,  n)  contains  data  for  boimdary  edges  m. 
ibdryCm,  1)  must  be  m. 

ibdry  (m,  2)  is  the  number  of  the  tri£Uigle  containing  the  boimdary  edge,  in  the  range 
1 to  nt. 

ibdry  (m,  3)  is  the  edge  number  of  the  boimdary  edge  within  the  triangle.  It  must 
be  1,  2,  or  3;  the  first  edge  is  opposite  the  first-named  vertex  in  the  itnode  array  for 
the  triangle. 

ibdry  (m,  4)  is  0 if  the  edge  is  straight,  and  is  the  curved  edge  number,  in  the  range 
1 to  nc,  if  the  edge  is  curved. 

ibdryCm,  n+4)  is  the  boundary  condition  type  for  the  n-th  PDE  on  this  boundary 
edge;  it  should  be  1 for  Dirichlet-type  boundary  conditions  amd  2 for  Neumann-type 
boundairy  conditions. 

3.4.9  USRIN3  is  the  third  input  routine.  It  is  called  to  obtain  parameters  for  the  solution 
process.  It  is  called  only  if  oldsol  is  .FALSE.  The  default  version  calls  I0GET(4)  to  get 
a file  number,  and  calls  FREEIN  to  read  the  values  from  the  file.  All  the  airguments  are 
outputs. 

subroutine  USRIN3(mauan,  maxr,  ittype,  Ivlref , Ivlusr,  numlvl, 

* ithrsh,  iprsw,  iwork) 

naxm  is  the  “m”  parameter  for  multi-level  iterations  [BaLnk79|. 
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mayr  is  the  “r”  parameter  for  multi-levei  iterations  [BankTO]. 
ittype  is  the  type  of  smoothing  iteration. 

0:  point  Gauss-Seidel. 

1:  conjugate  gradient. 

2:  block  Gauss-Seidel. 

Ivlref  is  the  maximum  pre-refinement  level  assigned  to  any  vertex. 

Ivlusr  is  the  maximiim  pre-refinement  level  used  in  USRTST. 
numlvl  is  the  number  of  mesh  levels  allowed, 
ithrsh  is  the  threshold  percentage  for  mesh  adapting, 
ithrsh  = 0 : do  uniform  refinement. 

ithrsh  > 0 : refine  triangles  whose  estimated  error  is  at  least  this  percent  of  the 
largest  estimated  error. 

ithrsh  < 0 : refine  triangles  whose  estimated  error  puts  them  in  the  worst 
abs  (ithrsh)  percent  of  all  triangles. 

iprsw  determines  how  much  printing  is  done  as  the  solution  is  being  found. 

0 : as  little  as  possible. 

1 : more. 

2 : reasonable  for  interactive  use. 

3 : far  too  much. 

iwork  is  the  number  of  words  of  stack  spew:e  to  use  for  solving.  B2DE  will  use  the 
smaJlest  of  iwork,  the  value  initialized  by  USRMAI,  and  the  actual  available  space. 

3.4.10  USRREG  should  return  the  regularization  parameter  for  the  n-th  PDE.  In  calculating 
the  Jacobian  matrix,  B2DE  uses  a„  + USRREG  (n)  to  allow  for  regularizing  nezirly-singular 
Jacobian  matrices.  The  default  version  returns  zero. 

real  function  USRREG (n) 
integer  n 

common  / us r com  / par  am 
real  param 

3.4.11  USRSAV  is  called  when  a solution  is  being  saved,  for  the  user  to  save  necessary 
problem  paurameters  not  known  to  B2DE.  The  default  version  does  nothing.  USRGET  should 
read  values  in  the  same  order  and  with  the  same  formats  as  USRSAV  writes.  USRSAV  is  always 
called  with  isout  equal  to  io(8). 

subroutine  USRSAV(isout) 
integer  isout 

isout  is  the  file  number  to  write  to. 

3.4.12  USRSHF  is  called  to  “shift,”  or  transform  PDE  variables  for  plotting.  All  the  argu- 
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ments  are  inputs  except  unew,  which  is  an  output.  The  default  version  does  nothing. 

subroutine  USRSHFCixtype , npde,  nv,  u,  unew,  vx,  vy) 

integer  ixtype,  npde,  nv 

real  u(npde,nv),  unew(nv),  vx(nv) , vy(nv) 

ixtype  is  the  index  of  the  transform  type. 

npde  is  the  number  of  PDEs. 

nv  is  the  number  of  vertices  in  the  current  finest  mesh. 

u is  the  solution  on  the  finest  mesh. 

vx  and  vy  are  the  coordinates  of  the  vertices. 

unew  is  the  transformed  solution.  It  may  be  any  desired  function  of  the  PDE  variables, 
the  coordinates,  or  anything  else. 

3.4.13  USRSTP  provides  the  user  an  opportunity  to  avoid  problems  such  as  overfiow  or  the 
square  root  of  a negative  number,  by  restricting  the  maximum  fraction  of  a Newton  step 
that  can  be  used.  B2DE  will  try  to  evaluate  residuals  at  (uold  -f-  t uatep),  where  t is  the 
value  returned  by  USRSTP.  USRSTP  should  return  a number  greater  than  0 and  less  than  or 
equal  to  1.  All  the  arguments  are  inputs.  The  default  version  returns  1. 

real  function  USRSTP (npde,  nv,  vx,  vy,  uold,  ustep) 
integer  npde.  nv 

real  vx(nv) , vy(nv) , uold (npde,  nv) , ustep(npde,  nv) 
npde  is  the  number  of  PDEs. 

nv  is  the  number  of  vertices  in  the  current  finest  mesh, 
vx  and  vy  are  the  coordinates  of  the  vertices, 
uold  is  the  ciirrent  solution, 
ustep  is  the  tentative  Newton  step. 

3.4.14  USRTST  allows  pre-refinement  of  the  initial  mesh  in  addition  to  that  provided  by 
the  array  Iref  in  USRIN2.  After  the  pre-refinement  scheduled  by  Iref , USRTST  is  called 
for  each  triangle,  and  should  return  .TRUE,  if  the  triangle  is  to  be  refined  at  least  one 
more  time.  All  triangles  are  tested  until  no  more  triangles  are  refined.  All  the  arguments 
are  inputs.  The  default  version  always  returns  . FALSE . 

logical  function  USRTST (x,  y,  Ivltri,  level 1) 
real  x(3) , y(3) 
integer  Ivltri , levell 

X and  y are  the  coordinates  of  vertices  of  the  triangle. 

Ivltri  is  the  level  of  refinement  of  this  triangle  so  far. 

levell  is  the  maximum  level  of  pre-refinement  allowed,  as  specified  by  USRIN3. 

3.4.15  USRWTS  is  called  when  the  mesh  is  being  refined;  it  should  return  the  relative  weight 
for  the  n-th  PDE.  If  any  of  the  weights  is  negative  euid  the  run  is  interewrtive,  the  driver 
asks  what  the  weights  should  be,  ignoring  the  values  returned  by  USRWTS.  The  default 
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version  always  returns  1. 

real  function  USRWTS(n) 
integer  n 

n is  the  PDE  number. 

4.  Writing  Your  Own  Driver  for  B2DE 

The  routines  callable  by  the  user  axe  described  in  several  groups:  DRVxxx,  SETxxx,  SOLxxx, 
PLTxxx,  PRTxxx,  TRMxxx,  USRxxx,  and  miscellaneously-named  routines. 

DRVxxx  routines  are  driver  programs,  and  can  be  modified  or  replaced  by  the  user.  DRVALL 
and  DRVSOL  are  described  in  Section  3.2,  and  DRVPLT  is  described  in  Section  3.3. 

DRVSEE  can  be  used  to  inspect  the  solution 

subroutine  DRVSEE (ausr) 
external  ausr 

DRVSEE  provides  the  menu  discussed  in  Section  3.3.3: 

Imprint , 2=plot , 3=continue 

SOLxxx  routines  do  the  solution  of  the  PDEs.  They  are  described  in  Section  4.2. 

PLTxxx  routines  do  plotting.  They  are  described  in  Section  4.3. 

PRTxxx  routines  do  printing.  They  are  described  in  Section  3.2.1. 

SETxxx  routines  set  parameters  used  for  plotting.  They  are  described  in  Section  4.1. 

TRMxxx  routines  are  terminal-dependent  programs  having  to  do  with  plotting.  They  axe 
described  in  Section  5.2. 

USRxxx  routines  provide  the  user  with  opportunities  to  do  things  without  requiring  detailed 
knowledge  of  the  interior  of  B2DE.  They  axe  described  in  Section  3.4. 

4.1  SETxxx  Routines 

4.1.1  SETCOL  sets  the  color  or  line  type  for  plotting.  It  saves  the  input  value,  then  calls 
the  terminal-dependent  routine  TRMCOL;  TRMCOL  is  described  in  Section  5.2.1. 

subroutine  SETCOL (newcol) 
integer  newcol 

newcol  is  the  new  color  or  line  type. 

4.1.2  SETERA  either  erases  the  screen  or  sets  the  fiag  that  determines  whether  the  screen  is 
erased  before  a new  plot  is  drawn.  On  some  terminals,  it  may  not  be  possible  or  desirable 
to  plot  without  erasing  the  screen. 

subroutine  seteraCif lag) 
integer  iflag 

If  iflag  is  1,  the  screen  will  not  be  erased  before  each  plot. 

If  iflag  is  2,  the  screen  will  be  erased  before  each  plot. 

If  iflag  is  3,  the  screen  is  erased  immediately. 
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Other  values  for  if  lag  do  nothing. 

4.1.3  SETMAG  sets  the  magnification  emd  the  x-  and  y-centers  for  the  magnified  plots. 
Magnification  applies  to  triangle  plots,  contour  plots,  and  fiow  line  plots.  It  also  applies 
to  printing  of  PDE  values  on  a grid. 

subroutine  setmagCxymag,  xctr,  yctr) 
real  xymag,  xctr,  yctr 

xymag  is  the  magnification.  If  xymag  is  less  than  1,  B2DE  uses  a magnification  of  1. 

xctr  and  yctr  are  the  coordinates  of  the  center  of  the  image  to  be  magnified,  in 
percentages  of  the  full  window  size.  For  example,  if  xctr  and  yctr  axe  each  50,  the 
center  of  the  image  will  be  magnified. 

4.1.4  SETWDW  sets  the  window  in  which  subsequent  plots  will  be  drawn. 

subroutine  setwdwCxmin,  xmax,  ymin,  ymax) 
real  xmin,  xmax.  ymin.  ymax 

xmin  and  xmax  aure  the  x-coordinates  of  the  window,  in  percentages  of  the  maxim\im 
window  size  as  given  by  TRMSCL.  Similarly,  ymin  and  ymax  axe  the  y-coordinates  of 
the  window.  For  exatmple,  if  xmin  amd  ymin  are  eau:h  50,  and  xmax  aind  ymax  are  each 
100,  the  plot  window  will  be  the  upper  right  quarter  of  the  plotting  area. 

4.1.  SETFPL  sets  a fiag  for  plotting  and  printing. 

subroutine  SETFPL (flag) 
logical  flag 

If  flag  is  .TRUE.,  subsequent  printing  and  plotting  will  be  done  to  a file. 

If  flag  is  .FALSE.,  subsequent  printing  and  plotting  will  be  done  to  the  terminal. 

4.1.  SETSHF  is  used  to  “shift”,  or  transform,  the  first  PDE  variable  for  inspecting  a 
user-defined  function  of  the  PDE  variables,  and  to  shift  it  back. 

logical  function  SETSHF (if lag) 
integer  if lag 

If  if  lag  is  not  0,  SETSHF  calls  USRSHF  to  get  new  values  to  put  in  place  of  the  values 
of  the  first  PDE,  exchzmges  the  values,  and  returns  for  plotting  or  printing,  if  lag  is 
passed  to  USRSHF  for  use  as  an  index. 

If  if  lag  is  0,  SETSHF  restores  the  values  saved  for  the  first  PDE. 

4.2  SOLxxx  Routines 

SOLxxx  routines  are  called  by  driver  programs  to  solve  the  PDEs.  They  are  described 
individually  below. 

4.2.1  SOLCHG  is  used  to  change  parameters  during  a solution.  There  are  three  types  of 
parameters.  First,  USRCHG  is  called  for  changing  user-supplied  pexameters.  Second,  if  there 
are  two  or  more  PDEs  in  the  problem,  one  or  more  PDE  vEiriables  can  be  held  fixed  during 
the  nonlinear  iterations.  Third,  B2DE  parixneters  maxm,  maxr,  ittype,  ithrsh,  and  iprsw 
can  be  changed. 
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subroutine  SOLCHG 

4.2.2  SOLCMP  is  used  to  compress,  or  merge,  the  top  two  levels  of  mesh  into  a single  mesh. 
Lower  level  meshes  are  not  affected.  SOLCMP  is  useful  if  a refinement  did  not  add  many 
new  vertices  to  the  mesh.  There  are  no  inputs  or  outputs. 

subroutine  SOLCMP 

4.2.3  SOLDIF  is  used  to  inspect  the  “difference”  in  a solution;  this  is  the  difference  between 
the  starting  solution  values  (initial  guess,  read-in  values  from  old  solution,  or  interpolated 
solution  from  next  coarser  mesh)  on  the  current  finest  mesh  ajid  the  current  solution 
values.  Norms  of  the  changes  are  printed,  and  DRVSEE  is  called  for  possible  printing  or 
plotting.  SOLDIF  is  useful  for  inspecting  the  changes  in  the  solution.  The  only  input  is  the 
a function;  there  are  no  outputs. 

subroutine  SOLDIF (ausr) 
external  ausr 

4.2.4  SOLINP  handles  most  of  the  details  of  input  for  a problem.  All  the  arguments  are 
outputs.  Their  meanings  were  discussed  in  Section  3.4.7,  under  USRINl. 

subroutine  SOLINP (oldsol,  nt.  nv.  nc . nb,  npde) 
logical  oldsol 
integer  nt . nv . nc . nb 

4.2.5  SOLITR  handles  the  nonlinear  iterations;  it  is  the  workhorse  routine.  The  inputs 
are  the  names  of  the  three  subroutines  describing  the  PDEs.  SOLITR  asks  for  values  of 
iteration  parameters  and  then  proceeds. 

subroutine  SOLITRCausr,  fusr,  gusr) 
external  ausr.  fusr,  gusr 

ausr,  fusr,  and  gusr  are  the  names  of  subroutines  that  describe  the  PDEs  and  the 
boundary  conditions. 

4.2.6  SOLREF  is  used  either  to  start  a new  solution  or  to  refine  the  mesh  for  an  existing 
solution. 

subroutine  SOLREF (ausr,  fusr.  gusr.  nevlvl) 
external  ausr.  fusr,  gusr 

ausr,  fusr,  and  gusr  are  the  names  of  subroutines  that  describe  the  PDEs  and  the 
boundziry  conditions. 

newlvl  should  be  0 to  get  started  on  the  same  level,  or  1 to  refine  to  the  next  higher 
level. 

4.2.7  SOLSAV  saves  the  solution  values  eind  the  data  structure  by  writing  them  to  a file 
in  a format  suitable  for  later  use  by  B2DE.  SOLSAV  calls  USRSAV  in  case  the  user  wants  to 
write  anything  for  later  retrieval  by  USRGET. 

subroutine  SOLSAV 

4.3  PLTxx  Routines 
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PLTxxx  routines  are  callable  by  the  user,  but  should  not  be  changed.  In  the  default  driver, 
only  the  corresponding  DRVxxx  routines  call  the  PLTxxx  routines.  The  five  routines  are 
PLTCON,  PLTFLW,  PLTPRO,  PLTSUR,  and  PLTTRI.  They  are  described  separtately. 

4.3.1  PLTCON  does  a contour  plot  of  a PDE  variable.  The  contour  values  are  printed  to 
the  standard  output,  but  do  not  appear  on  the  plot.  The  plot  uses  the  current  color  (or 
color  palette  for  filled  contours) , window,  and  magnification. 

subroutine  PLTCON (nth,  ntype,  fill,  istcol,  nlev,  ulev) 
integer  nth,  ntype,  istcol,  nlev 
real  ulev(*) 
logical  fill 

nth:  Plot  contours  of  the  n-th  PDE. 

ntype:  If  ntype  is  1,  plot  contours  of  u„(x,  y);  otherwise  plot  contours  of  | Vun(i,  y)  |. 

fill:  If  fill  is  .TRUE. , the  spaces  between  contours  are  filled  in;  otherwise  only  the 
contour  lines  lire  drawn. 

istcol  is  the  starting  color  for  filling  spaces  between  contours,  istcol  is  ignored  if 
fill  is  .FALSE. 

nlev:  If  nlev  is  positive,  plot  nlev  contours,  using  automatic  scaling,  spacing  the 
contours  evenly  between,  but  not  including,  the  maximum  and  minimum  values  of 
the  plot  variable;  ulev  is  ignored  by  PLTCON.  If  nlev  is  negative,  plot  -nlev  contours, 
using  the  contour  values  in  the  ulev  array.  If  nlev  is  zero,  just  find  and  print  the 
minimum  and  maximum  values  of  the  plot  variable. 

4.3.2  PLTFLW  draws  flow  lines  for  a pde  variable.  The  amount  of  positive  and  negative 
flow  is  printed  to  the  standard  output.  The  plot  uses  the  current  colors,  window,  and 
magnification. 

subroutine  PLTFLW(axy,  nth,  npf , nflow,  xyval,  rp) 
real  xyval (♦) 
external  axy 

axy  is  the  name  of  the  a subroutine, 
nth  is  as  described  under  PLTCON. 
npf  is  the  number  of  points  defining  the  crossing  line, 
xyval  is  an  airray  of  the  points,  first  the  npf  z- values,  then 
nflow  is  the  number  of  flow  lines  to  be  plotted. 

4.3.3  PLTPRO  draws  profile  plots,  the  outline  of  slices  through  a 
uses  the  current  color  emd  window,  but  uses  magnification  of  1. 

subroutine  PLTPRO (nth,  nprof,  dir) 
integer  nth,  nprof 
real  dir (2) 

nth  is  as  described  under  PLTCON. 


the  npf  y- values. 


projection  plot.  The  plot 
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nprof  is  the  number  of  profiles. 

dir  are  the  indices  of  the  slicing  plane,  as  discussed  in  Section  3.3.4. 

4.3.4  PLTSUR  draws  a surface  projection  plot  of  a PDE  variable.  The  plot  uses  the  current 
color  and  window,  but  uses  magnification  of  1. 

subroutine  PLTSURCnsurf , level,  nth,  ntype,  vievrpt,  ngrid, 

* nlev,  ulev) 

integer  nsurf,  level,  nth,  ntype,  ngrid (2) , nlev 
real  viewpt(3),  ulev(*) 

nsurf  controls  what  is  plotted  on  the  surface.  If  nsurf  is  1,  trizingles  of  the  specified 
level  are  plotted.  If  nsiirf  is  2,  contours  are  plotted.  If  nsurf  is  3,  x-y  grids  are 
plotted. 

level  is  the  level  of  solution  to  be  plotted, 
nth  and  ntype  a^e  as  described  under  PLTCON. 

viewpt  contains  the  (r,  y,  2)  components  of  the  view  point.  The  surface  projection 
will  appear  as  if  it  were  at  the  origin,  seen  from  this  point. 

nlev  amd  ulev  are  used  if  only  nsurf  is  2.  Their  meauings  are  as  given  under  PLTCON. 

4.3.5  PLTTRI  draws  a plot  of  one  triangulation.  The  plot  uses  the  current  color,  window, 
and  magnification. 

subroutine  PLTTRI (level) 
integer  level 

level  is  the  level  of  the  mesh  to  be  plotted. 

4.4  Miscellaneous  routines 

4.4.1  lOGET  is  used  to  get  the  vadue  of  one  of  the  input/output  pau’ameters.  See  the  USRMAI 
section  for  a list  of  the  meatnings. 

integer  function  ioget (ionum) 
integer  ionum 

ionum  is  the  index  of  the  paxauneter  whose  value  is  desired. 

4.4.2  lOSET  is  used  to  set  am  input/output  pairameter.  See  the  USRMAI  section  for  a list 
of  the  meanings. 

subroutine  ioset (ionum,  ioval) 
integer  ionum.  ioval 

ionum  is  the  index  of  the  pairauneter  to  be  set. 

ioval  is  the  value  to  set. 

4.4.3  FREEIN  is  a subroutine  to  do  free-format  input.  The  calling  prograim  asks  for  some 
real  numbers  amd/or  integers  to  be  read  from  am  input  file.  FREEIN  reads  a line  at  a time, 
parses  it,  evaluates  the  numbers,  and  does  error  checking.  If  an  error  occurs  and  the  session 
is  interactive,  FREEIN  prompts  the  user  for  a new  line.  Items  are  sepau-ated  by  spaces;  lines 
aue  read  until  all  the  requested  items  are  obtained.  If  there  are  more  items  on  a line  than 
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the  call  to  FREEIN  requires,  the  excess  is  ignored.  A * or  # character  majks  the  beginning 
of  a comment;  the  character  and  the  rest  of  the  line  are  ignored. 

subroutine  FREEIN (ioread,  ioerr,  nitems,  itype,  ibuf f , rbuff) 
integer  ioread,  ioerr,  nitems,  itype (nitems) , ibuf f (nitems) 
real  rbuff (nitems) 

ioread  is  the  file  number  from  which  to  read  the  data, 
ioerr  is  the  file  number  on  which  to  write  error  messages, 
nitems  is  the  number  of  items  to  be  read. 

itype(j)  is  the  type  of  the  j-th  item,  1 for  an  integer  and  2 for  a real  number. 

ibuf f and  rbuff  are  arrays  into  which  the  data  is  put.  If  itype  (j ) is  1,  the  j-th  item 
is  interpreted  as  an  integer  and  put  into  ibuff  (j);  if  itype(j)  is  2,  the  j-th  item  is 
interpreted  as  a real  number  and  put  into  rbuff  ( j ) . 

4.4.4  FLTGET  is  provided  for  convenience;  it  returns  one  real  value  typed  at  the  terminal, 
using  FREEIN. 

real  function  fltget(junk) 
j unk  is  a dummy  parameter. 

4.4.5  INTGET  is  provided  for  convenience;  it  returns  one  integer  value  typed  at  the  terminal, 
using  FREEIN. 

integer  function  intget(junk) 
junk  is  a dummy  parameter. 

4.4.6  EVAL  evaluates  PDE  variables  for  passing  to  other  programs. 

subroutine  EVAL (kind.  nth.  npde,  npts,  x.  y.  u.  dudx,  dudy) 
real  x(npts) , y(npts) , u(npts) , dudx(npts) , dudy(npts) 

If  kind  is  1,  only  the  PDE  variable  is  evaluated  and  returned  in  u;  dudx  zmd  dudy  are 
not  referenced.  If  kind  is  2,  only  the  gradient  of  the  PDE  variable  is  evaluated  and 
retiirned  in  dudx  and  dudy;  u is  not  referenced.  If  kind  is  3,  both  the  PDE  variable 
and  its  gradient  are  evaluated. 

The  nth  PDE  variable  is  evaluated. 

npde  is  the  number  of  PDE  variables. 

npts  is  the  number  of  points  at  which  to  evaluted  the  nth  PDE. 

X and  y are  the  points  at  which  to  evaluate  the  nth  PDE. 

If  kind  is  1 or  3,  the  PDE  values  are  returned  in  u. 

If  kind  is  2 or  3,  the  PDE  gradient  values  zu'e  returned  in  dudx  and  dudy. 

5.  Installing  B2DE  on  Your  Computer 

Most  of  B2DE  is  standard  Fortran  77;  it  should  not  need  changes  by  the  user  euid  should 
not  be  changed.  All  known  machine  dependencies  have  been  isolated  for  easy  changing. 
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5.1  The  PORT  Stack 


Working  space  is  managed  inside  B2DE  using  the  PORT  stack  [Fox78].  The  size  of  the 
stack  determines  how  large  a problem  can  be  solved.  The  stack  takes  advantage  of  a 
feature  which  is  in  most  Fortrans,  but  which  is  not  in  the  Fortran  77  standard:  the  size  of 
a labelled  common  block  is  determined  by  the  size  in  the  main  program,  or  by  the  size  in 
the  first  program  loaded.  To  change  the  size  of  the  stack,  only  the  main  program  needs  to 
be  edited  and  re-compiled.  The  following  block  of  code  occurs  many  times  within  B2DE, 
defining  the  stack  common,  cstak,  to  be  a nominal  size. 

parameter  (irsize=1000) 

parameter  (idsize-  500) 

common/cstak/dstak 

double  precision  dstak(idsize) 

real  r(irsize) 

integer  i(irsize) 

logical  IstedcCirsize) 

equivalence  (dstak(l) .r(l)) . (dstak(l) .i(D) , (dstak(l) .Istak(l) ) 
save  /cstedc/ 

Your  main  program  will  define  the  true  size  of  the  stack. 

5.2  Machine-Dependent  Programs 

To  be  as  portable  as  possible,  B2DE  uses  the  machine  constants  from  the  PORT  library 
[Fox78j.  The  machine  constant  functions  IlMACH  and  RIMACH  must  be  changed  as  appropri- 
ate for  your  computer.  The  distributed  version  is  correct  for  the  DEC  VAX  11/700  series 
and  the  VAX/VMS  operating  system.  Constants  for  many  other  computers  are  in  the 
routines,  but  commented  out.  See  [Fox78]  for  instructions  on  changing  machine  constants 
for  computers  not  listed. 

B2DE  \ises  two  timing  routines.  Each  should  calculate  an  elapsed  amoimt  of  time  used  in 
milliseconds.  ISTIME  is  called  at  the  beginning  of  an  activity  to  be  timed,  and  IFTIME  is 
called  at  the  conclusion.  The  difference  is  used  as  the  elapsed  time. 

integer  function  iftime  ( j\ink  ) 

junk  is  a dummy  parameter. 

integer  f\mction  istime  ( junk  ) 

junk  is  a dummy  parameter. 

5.3  Terminal-Dependent  Programs 

TRMxxx  routines  are  terminal-dependent.  The  default  versions  work  for  DEC  VT125  ter- 
minals iind  Hewlett  Packard  264x  and  262x  terminals.  New  versions  must  be  written  for 
other  terminals.  Some  TRMxxx  routines  are  not  called  if  plotting  to  a file  is  being  done,  and 
therefore  are  not  necessary  if  only  batch  operation  is  planned,  although  dummy  versions 
must  be  provided  to  satisfy  most  operating  systems. 

5.2.1  TRMCOL  should  set  the  color  or  the  line  type  for  the  plots.  If  the  terminal  does  not 
support  color  or  alternate  line  types,  TRMCOL  need  not  do  anything.  TRMCOL  is  called  by 
SETCOL,  which  is  called  by  DRVCOL;  DRVCOL  only  allows  non-negative  values  for  the  color. 


31 


subroutine  TRMCOL(newcol) 
integer  newcol 

newcol  is  the  new  color  or  line  type. 

5.2.2  TRMDRW  should  draw  a line  from  the  last  point  reached  by  TRMDRW  or  TRMMOV  to  the 
point  (ix,  iy).  The  point  (ix,  iy)  is  guaranteed  to  be  within  the  limits  set  by  TRMSCL. 
TRMDRW  is  not  called  if  plotting  to  a file  is  being  done.  The  line  should  be  drawn  in  the 
current  color  or  line  type  as  last  set  by  TRMCOL. 

subroutine  TRMDRW (ix,  iy) 
integer  ix,  iy 

ix  and  iy  are  the  screen  coordinates  of  the  new  point. 

5.2.3  TRMMOV  should  move  the  current  point  to  the  point  (ix,  iy).  No  line  should  be 
drawn.  The  point  (ix,  iy)  is  guaranteed  to  be  within  the  limits  set  by  TRMSCL.  TRMMOV 
is  not  called  if  plotting  to  a file  is  being  done. 

subroutine  TRMMOV (ix,  iy) 
integer  ix,  iy 

ix  and  iy  eire  the  screen  coordinates  of  the  new  point. 

5.2.4  TRMPAL  changes  the  color  palette  available  to  the  user.  It  has  no  arguments, 
subroutine  TRMPAL 

5.2.5  TRMPLT  is  a plotting  utility.  Its  function  depends  upon  its  argument,  as  described 
below. 

subroutine  TRMPLT (if lag) 
integer  if lag 

if  lag  = 1:  TRMPLT  should  do  whatever  is  necessary  to  initialize  all  plotting  for  the 
session. 

if  lag  = 2:  TRMPLT  should  do  whatever  is  necessary  to  begin  a plot,  amd  should  erase 
the  screen. 

if  lag  = 3:  TRMPLT  should  do  whatever  is  necessary  to  begin  a plot,  but  should  not 
erase  the  screen.  This  is  not  possible  on  some  terminals. 

if  lag  = 4:  TRMPLT  should  do  whatever  is  necessary  to  finish  a plot.  If  the  session  is 
interactive,  and  if  plots  are  done  on  the  same  screen  as  text,  TRMPLT  should  wait  for 
a signal  from  the  user  before  returning. 

iflag  = 5:  TRMPLT  should  do  whatever  is  necessary  to  finish  aill  plotting  for  the 
session. 

5.2.6  TRMSCL  should  give  the  minimum  and  maximum  screen  coordinates  for  terminal  type 
iterm.  If  iterm  is  zero  or  am  unknown  type,  the  values  for  plotting  to  a file  should  be 
returned. 

integer  function  TRMSCL (iterm.  iflag) 
integer  iterm.  iflag 

iflag  = 1:  TRMSCL  should  return  the  minimum  x value. 


32 


if  lag  = 2:  TRMSCL  should  return  the  maximum  x value, 
if  lag  = 3:  TRMSCL  should  return  the  minimum  y value, 
if  lag  = 4:  TRMSCL  should  return  the  maximum  y value. 

6.  Examples 

6.1  Laplace’s  Equation 

The  equation  is  Laplace’s  equation,  the  simplest  elliptic  PDE.  The  region  is  a square  with 
two  comers  roimded  off,  as  in  Fig.  1. 

V • Vu  = 0 


The  boimdary  conditions  2ire  u = 1 on  edges  1 zmd  2,  u = 0 on  edge  5,  zmd  duj dv  — Q on 
sides  3,  4,  and  6.  The  solution  is  singular  at  vertices  1 and  6. 

Appendix  A is  a program  to  solve  this  example;  Appendix  B is  an  input  data  file  suitable 
for  use  with  the  default  driver;  and  Appendix  C is  a sample  run. 

A typical  initiad  triangulation  supplied  by  a user  is  shown  in  Figure  2.  Sides  1 and  4 are 
defined  as  circular  arcs,  but  the  curvature  does  not  show  up  imtil  the  sides  are  refined. 
Figure  3a  is  the  actual  coarsest  triangulation  used;  the  input  data  file  specifies  initial 
refinement  arotmd  the  singularities.  Figures  3b  amd  3c  are  refined  meshes  generated  by 
B2DE  in  the  sample  run.  Note  that  the  adaptive  refinement  concentrated  new  triangles 
neax  the  singularities.  Figure  4 shows  other  plots  froni  the  sample  run:  4a  is  a contour 
plot;  4b  is  a fiow  line  plot;  and  4c  is  a profile  plot. 

6.2  A Nonlinear  Example 

This  is  a concocted  set  of  three  nonlinear  elliptic  PDEs.  The  solutions  are  smooth,  and  the 
adaptive  mesh  generates  uniform  triangulations.  The  region  is  the  rectangle  — 1 < a:  < 1, 
0 < y < 1. 


V • (V U2)  = 0 
V • (U1U3VU3)  = 0 

The  boundary  conditions  are  all  Dirichlet: 


U1U3U3  = 1 


U3  = e* 
U3  = 
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The  solution  is 


U2  = e® 

Us  = 

Appendix  D is  a program  to  solve  this  example;  Appendix  E is  an  input  data  file  suitable 
for  use  with  the  default  driver;  and  Appendix  F is  a sample  run. 

A typical  initial  triangulation  supplied  by  a user  is  shown  in  Figure  5a.  Figure  5b  is  the 
actual  coarsest  triangulation  used;  the  input  data  file  specifies  initial  refinement  around  the 
left  comer.  (There  is  no  singularity  in  the  problem;  this  initial  refinement  is  unnecessary.) 
Figures  5c  zuid  5d  are  refined  meshes  generated  by  B2DE  in  the  sample  run.  Note  that  the 
adaptive  refinement  made  the  mesh  uniform;  the  incorrect  initial  mesh  has  been  overcome. 

Figure  6 shows  other  plots  from  the  sample  nm:  6a  is  a surface  plot  showing  triangles;  6b 
is  a surfax:e  plot  showing  contours,  with  hidden  lines  not  removed;  and  6c  is  a surface  plot 
showing  grid  lines. 
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1.  Region  for  Laplace’s  equation  example.  The  boundary  conditions  are  u = 1 on  edges  1 
and  2,  u = 0 on  edge  5,  and  duifdv  = 0 on  sides  3,  4,  and  6. 

2.  User-supplied  triangulation  of  region.  The  curvature  of  sides  1 and  4 will  show  up  in 
refinements  of  the  region. 

3a.  Initial  triangulation  of  region.  Vertices  1 and  6 were  specified  as  refinement  level  3; 
other  vertices  were  level  1. 

3b.  Triangulation  after  first  adaptive  refinement. 

3c.  Triangulation  after  second  adaptive  refinement. 

4a.  Contours  of  the  approximate  solution  on  the  mesh  of  Fig.  3c.  The  contours  ajre  imiformly 
spaced. 

4b.  Flow  line  plot  of  the  approximate  solution  on  the  mesh  of  Fig.  3c.  Equal  fiux  fiows 
between  each  pair  of  fiow  lines. 

4c.  Profile  plot  of  the  approximate  solution  on  the  mesh  of  Fig.  3c,  with  slicing  direction 

1 0. 

5a.  User-supplied  triangulation  of  region  for  example  2. 

5b.  Initial  triangulation  of  region.  Vertex  4 was  specified  as  refinement  level  3;  other  vertices 
were  level  1. 

5c.  Triangulation  edter  first  adaptive  refinement. 

5d.  Triangulation  after  second  adaptive  refinement. 

6a.  Surface  plot  of  error  in  PDE  2,  on  mesh  of  Fig.  5d.  Triangles  are  shown,  ajid  hidden 
lines  are  removed. 

6b.  Surface  plot  of  error  in  PDE  2,  on  mesh  of  Fig.  5d.  Contours  are  shown,  and  hidden 
lines  are  not  removed. 

6c.  Surfax:e  plot  of  error  in  PDE  2,  on  mesh  of  Fig.  5d.  Grid  lines  are  shown,  21  along  x 
and  11  along  y,  and  hidden  lines  are  removed. 
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Appendix  A.  Program  to  Solve  Laplace’s  Equation  Example 
c 

c main  program  to  drive  nonlinear  multi -variable  code  to  solve 
c Laplace’s  equation,  adapted  from  umain.for 
c 

c the  user-supplied  routines  that  determine  the  PDEs 
c 

external  alap,  flap,  glap 
c 

c i/o  value  array 
c 

integer  io(9) 
c 

c the  stack  must  be  declared  and  initialized  in  the  main  program, 
c it  is  ’save’ed  for  safety’s  sed^e.  the  equivalence  is  required  on 
c some  con^uters  to  get  the  word  alignment  correct, 
c 

parameter  (irsize  = 100000) 
parameter  (idsize  = 60000) 

common/cstedc/dstedc 
double  precision  dstak(idsize) 
real  rstak(irsize) 
equivalence  (dstedcCl)  . rstak(l)) 
save  /cstak/ 
c 

c initialize  the  stack  to  tell  it  the  true  total  size 
c 

call  istkinCirsize,  3) 
c 

c specify  i/o  values 
c 

c io(l)  = file  number  for  reading  from  the  terminal 

c io(2)  » file  number  for  writing  to  the  terminal 

c io(3)  * terminal  type,  used  by  graphics  programs 
c 

io(l)  = 5 
io(2)  * 6 
io(3)  » 126 
c 

c io(4)  = file  number  for  reading  input  data; 
c only  used  by  usrinl . usrin2 , and  usrin3 

c io(6)  = file  number  for  writing  plot  data,  if  not  to  terminal 

c io(6)  * file  number  for  writing  print  data,  if  not  to  terminal 

c 

io(4)  = 20 
io(6)  * 21 
io(6)  = 22 
c 
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c io(7)  = file  number  for  reading  old  solution  data 

c io(8)  = file  number  for  writing  new  solution  data 

c io(9)  = file  number  for  writing  scratch  data 

c For  reading  aind  writing  solution  data,  the  data  is  written  as 

c unformatted,  or  binary,  if  the  io  value  is  positive,  and  as 
c formatted  if  the  io  value  is  negative, 
c 

io(7)  = -23 
io(8)  = -24 
io(9)  = 25 
c 

c install  the  i/o  values 
c 

do  10  j = 1.  9 

call  iosetCj,  io(j)) 

10  continue 
c 

c open  data  files, 
c 

openCunit  = io(4),  file  = ’inputl.dat’,  status  = ’old’) 
openCunit  = io(5) , file  = 'plot.daf  , status  = ’unknown’) 
openCunit  * io(6),  file  = ’print.dat’  , status  = ’unknown’) 
c 

if  ( io(7)  .gt.  0)  then 

openCunit  = ioC7),  file  » ’insol.dat’,  status  = ’unknown’, 

* form* ’unformatted’) 
else 

openCunit  = -ioC7) , file  = ’insol.dat’,  status  = ’unknown’, 

* form* ’formatted’) 
endif 

c 

if  C ioC8)  .gt.  0)  then 

openCunit  * ioC8) , file  * ’outsol.dat’,  status  = ’unknown’, 

* form* ’unformatted’) 
else 

openCunit  = -ioC8),  file  = ’outsol.dat’,  status  = ’unknown’, 

* form* ’formatted*) 
endif 

c 

openCunit  * ioC9) , file  * ’scrtch.dat’,  status  * ’\mknown’, 

* form* ’unformatted’) 
c 

c teike  it  away  . . . 
c 

call  drvallCalap,  flap,  glap) 
c 

c close  the  files,  just  to  be  safe 
c 
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close (io(4)) 
close(io(5)) 
close(io(6)) 
close (labs (io (7))) 
close(iabs(io(8) )) 
close (io (9)) 
c 

stop 

end 

c 

c 

subroutine  alap(npde,  ndim,  dopde,  x,  y,  i,  u,  dudx,  dudy, 

* a,  dadu,  dadux,  daduy) 
real  x,  y,  u(*) , dudx(*) , dudy(*) , 

* a(*) , daduCndim,*) , daduxCndim,*) , daduy(ndim, ♦) 
logical  dopde (npde) 

c 

c inputs: 

c npde  xthe  number  of  pdes  in  the  problem 

c ndim  the  first  dimension  of  arrays  dadu.  dadux,  eind  daduy 

c dopde  dopde (n)  is  .TRUE,  if  the  nth  is  active, 

c 7 the  coordinates  of  the  point  at  which  to  evaluate 

c i the  macro  triangle  number  that  the  point  is  in 

c u the  current  pde  values  at  the  point 

c dudx. dudy  the  gradient  of  the  pdes  at  the  point 

c 

c outputs : 

c a the  values  of  the  a coefficient 

c dadu  the  first  derivatives  of  a with  respect  to  u 

c dadux  the  first  derivatives  of  a with  respect  to  du/dx 

c daduy  the  first  derivatives  of  a with  respect  to  du/dy 

c 

a(l)  * 1. 
daduCl,  1)  * O.eO 
daduxCl,  1)  » O.eO 

daduy (1,  1)  =*  O.eO 

c 

return 

end 

c 

c 

subroutine  flap (npde,  ndim,  dopde,  x,  y,  i,  u,  dudx,  dudy, 

♦ f,  dfdu,  dfd\ix,  dfduy) 
real  x,  y,  u(*) , dudx(*) , dudy(*) , 

♦ f(*).  dfdu(ndim, *) , df dux (ndim, *) , dfduy (ndim, •) 
logical  dopde (npde) 

c 

c 
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c inputs: 

c npde  xthe  number  of  pdes  in  the  problem 

c ndim  the  first  dimension  of  airrays  dadu,  dadux,  and  daduy 

c dopde  dopde(n)  is  .TRUE,  if  the  nth  is  active, 

c 3C,  y the  coordinates  of  the  point  at  which  to  evaluate 

c i the  macro  triaoigle  number  that  the  point  is  in 

c u the  current  pde  values  at  the  point 

c dudx.dudy  the  gradient  of  the  pdes  at  the  point 

c 

c outputs: 

c f the  values  of  the  f coefficient 

c dfdu  the  first  derivatives  of  f with  respect  to  u 

c dfdux  the  first  derivatives  of  f with  respect  to  du/dx 

c dfduy  the  first  derivatives  of  f with  respect  to  du/dy 

c 

f(l)  = O.eO 
dfdu  (1.  1)  = O.eO 

dfduxCl,  1)  = O.eO 

dfduyCl,  1)  » O.eO 

c 

return 

end 

c 

subroutine  glapCnpde,  ndim,  dopde.  x.  y,  i,  j,  u.  g,  dgdu) 
real  x,  y,  u(*) , g(*) . dgduCndim,*) 
logical  dopde (npde) 
c 

c inputs: 
c npde 

c ndim 

c dopde 

c X.  y 

c i 

c j 

c u 

c 

c outputs : 

c g . 

c dgdu 

c 
c 

c square  with  2 round  comers 
c 

if  (i.eq.l  .or.  i.eq.2)  then 
g(l)  = u(l)  - 1.0 
dgdud,  1)  » 1.0 
else  if  (i  .eq.  5)  then 
g(l)  » u(l) 


xthe  number  of  pdes  in  the  problem 

the  first  dimension  of  arrays  dadu,  dadux,  and  daduy 
dopde (n)  is  .TRUE,  if  the  nth  is  active, 
the  coordinates  of  the  point  at  which  to  evaluate 
the  macro  triangle  number  that  the  point  is  in 
the  side  of  the  macro  triangle  that  the  point  is  on 
the  current  pde  values  at  the  point 

the  values  of  the  g coefficient 

the  first  derivatives  of  g with  respect  to  u 
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dgduCl.  1)  = 1.0 
else 

g(l)  = 0.0 
dgduCl,  1)  =0.0 
endif 

return 

end 
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Appendix  B.  File  “inputl.dat”  for  Laplace’s  Equation  Example 

# iziput.dat  for  Laplace  example 

It  square  with  two  rounded  comers 
# old  solution  flag 
7 2 6 1 # nt  nv  nc  nb  npde 

# vert  Ivl  X y 

1 3 -0.5  0.0 

2 1 0.0  0.5 

3 1 0.5  0.5 

4 1 0.0  0.0 

5 1 0.5  0.0 

6 3 0.0  -0.5 

7 1 -0.5  -0.5 

# curved  edges 

1 .35355339  -.35355339 

2 -.35355339  .35355339 

# triangles 
112  4 

2 2 3 4 

3 3 4 5 

4 4 5 6 

5 4 6 7 

6 14  7 . 

11321  # be  tri  edge  curv  be 

2 2 3 0 1 

3 3 2 0 2 

4 4 112 

5 5 10  1 

6 6 2 0 2 

2 8 0 # maxm  maxr  ittype 

4 17  It  Ivlref  Ivlusr  maxi 
-30  # ithresh 

2 # iprsw 

80000  # work  space 
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Appendix  C.  Transcript  of  interactive  session  for  Laplace’s  Equation  Example 
Output  from  the  computer  is  in  typewriter  font.  Numbers  typed  by  the  user  are  in 
boldface  and  have  a • at  the  beginning  of  the  line. 


refsol  entered  at  level  1 7 vertices 

grid  0.040  sec.  26  vertices 

asmbnl  0.030  sec. 

l=print,  2=plot,  3=iter,  4=refine,  5=comp,  6=save,  7=diff,  8=end,  9=change 

• 3 

enter  niter,  ieps.  relerr,  abserr,  jsmax,  paraun 

• 1 0 0 0 1 0 
Iter  1 

1 residual  L2  norm  4.397E-01 
factor  0.040  sec.  26  eqs. 

asmbnl  0 . 250  sec . 

mg  0.000  sec.  GE 

1 newton  norm, solution  norm, relerr  l.OOOE+00  O.OOOE+00  l.OOOE+00 

0=continue,  l=print,  2=plot 

• 0 

0 t 0.000000  L2(rhs)  4.397E-01  0.110  sec. 

1 t 1.000000  L2(rhs)  1.464E-07  0.110  sec. 

l»print,  2=plot,  3=iter,  4=refine,  5=comp,  6=save,  7=diff,  8=end,  9=change 

• 4 

refsol  re-entered  at  level  2 26  vertices 

asmbnl  0 . 030  sec . 

1 estimated  maximum  error  8.01E-03 
estimated  individual  PDE  errors 
0=continue , l=print , 2*plot 

• 0 


level  1 pde  1 \morm,  enorm  7.06E-01 
min  error  / max  error  1.44  percent 

4 of  30  triangles  ( 13.33)  above  50.00  percent 

9 of  30  triangles  ( 30.00) 

8 of  30  triamgles  ( 26.67) 

refine  0.300  sec.  level  2 

l=print,  2=plot,  3=iter,  4=refine, 

• 5 


2.86E-02  eat.  digits  1.39 


above  25.00  percent 
above  37.60  percent 
52  vertices 

5=comp,  6=save , 7=diff,  8»end,  9=change 


compressed  to  level  1 

l=print,  2=plot,  3=iter,  4=refine,  5=comp,  6=save,  7=dif f , 8=end,  9»change 

• 3 


enter  niter,  ieps,  relerr,  abserr,  jsmax.  param 
• 100010 
Iter  1 

1 residual  L2  norm  3.650E-02 
factor  0.080  sec.  52  eqs. 

asmbnl  0 . 480  sec . 

mg  0.000  sec.  GE 

1 newton  norm. solution  norm, relerr  7.664E-02  l.OOOE^OO  7.564E-02 
0=continue,  l=print,  2=plot 
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• 0 

0 t 0.000000  L2(rhs)  3.650E-02  0.260  sec. 

1 t 1.000000  L2(rhs)  7.992E-08  0.240  sec. 

l=print,  2=plot,  3=iter,  4=refine,  6=comp,  6=save,  7=diff,  8=end,  9=ch2mge 

• 4 


ref sol  re-entered  at  level  2 52  vertices 

asmbnl  0.050  sec. 

1 estimated  maximum  error  2.64E-03 
estimated  individual  PDE  errors 
0=continue,  l=print,  2=plot 
• 0 


level  1 pde  1 
min  error  / max 
9 of  72 

26  of  72 

14  of  72 

19  of  72 

23  of  72 

20  of  72 

21  of  72 

22  of  72 

refine  0.720 
l=print , 2=plot 
• 3 


unorm,  enorm  7.04E-01 
error  5.91  percent 
triangles  ( 12.50)  above  50.00  percent 


1.04E-02  est.  digits  1.83 


triangles  ( 36.11) 
triangles  ( 19.44) 
triangles  ( 26.39) 
triangles  ( 31.94) 
triangles  ( 27.78) 
trieoigles  ( 29.17) 
triangles  ( 30.56) 
sec . level  2 
, 3=iter,  4=refine, 


above  25.00  percent 
above  37.50  percent 
above  31.25  percent 
above  28.13  percent 
above  29.69  percent 
above  28.91  percent 
above  28.52  percent 
103  vertices 
5=comp,  6=save,  7=diff, 


8=end,  9=change 


enter  niter,  ieps,  relerr,  abserr,  jsmax,  param 
• 1 10  0 0 1 0 
Iter  1 

1 residual  L2  norm  1.761E-02 


factor 

0.080 

sec . 52 

eqs . 

asmbnl 

1.240 

sec . 

mg  2 

7.619E-01 

3.472E-02 

4.557E-02 

mg  2 

3.472E-02 

4.404E-03 

1.268E-01 

mg  2 

4.404E-03 

5.974E-04 

1.357E-01 

mg  0.160  sec.  GS  reduction  7.84E-04 

1 newton  norm, solution  norm, relerr  5.398E-02  l.OOOE+00  5.398E-02 
0=*continue,  l=print,  2*plot 


• 0 

0 t 0.000000  L2(rhs)  1.761E-02  0.510  sec. 

1 t 1.000000  L2(rhs)  1.659E-06  0.500  sec. 

l=print,  2*plot,  3=iter,  4*refine,  5=comp,  6=save,  7=diff,  8=end,  9=change 


• 1 


l=8ummary , 2»timing,  3=storage,  4=vertices,  5=pde  variables,  6=settings 

• 1 


n Ivl  vertices  triangles  unorm  enorm  digits 

11  52  114  7.04E-01  1.04E-02  1.83 

1 2 103  258 

l=summary,  2=timing,  3=storage,  4=vertices,  5=pde  variables,  6=settings 

• 2 
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timing  distribution 


plot 

seconds 

0.000 

assembly 

2.080 

factor 

0.200 

mg 

0.160 

rhs  check 

1.730 

refine 

1.060 

total 

5.230 

l=summary . 

2=timing , 

• 3 

103 

vertices 

2094 

vertices 

263 

trieuigles 

6583 

trieingles 

l=3ummary , 

2=timing , 

• 6 

percent 

0.00 

39.77 

3.82 

3.06 

33.08 

20.27 

100.00 

3=storage,  4=vertices, 

used. 

allowed. 

used. 

allowed. 

3=storage,  4=vertices, 


5=pde  variables , 6=setting3 


6=pde  variables,  6=setting3 


l=color,  2=erase,  3=magnify,  4=window,  5=file  print/plot,  6=transform 

• 3 


enter  magnification 

• 3 

enter  x-center  and  y-center  (0  to  100) 

• 50  0 

l=color,  2=erase,  3=magnify,  4=window,  6=file  print/plot,  6=transform 

• .0 

l«a\nranary,  2=timing,  3=storage,  4=vertices,  6=pde  variables,  6=“settings 

• 5 

l=on  grid.  2=at  vertices,  3=on  grid  as  scaled  integers 

• 3 

enter  pde  no.,  l=value  only,  2=gradient  only,  3=  both 

• 11 


nx,  ny 

• 77 

X limits  -1.66667E-01  1.66667E-01 

y limits  -6.00000E-01  -1.66667E-01 

pde  1 

scaled  values  from  -6.086E-08  to  6.492E-01  scale  5.492E-01 


90 

90 

91 

93 

96 

97 

100 

78 

78 

80 

82 

86 

89 

92 

64 

66 

67 

70 

76 

80 

86 

49 

60 

54 

68 

66 

71 

77 

34 

36 

39 

44 

64 

62 

70 

17 

19 

21 

28 

42 

66 

66 

0 

0 

0 

0 

9999 

9999 

9999 

l=on 

• 0 

grid. 

2=at 

vertices , 

3“on 

grid  as  scaled  integers 

l=summary , 

• 0 

2-timing, 

3=storage 

, 4=vertices,  6=pde  variables. 
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l=print,  2=plot,  3=iter,  4=refine,  5=comp,  6=save,  7=diff,  8=end,  9=change 

• 2 

I=tri2uigle8 , 2=conto\irs,  3=surface,  4=profiles,  6=flow  lines,  6=settings 

• 6 

l=color,  2=erase,  3=magnify,  4=window,  6=file  print/plot,  6=transform 

• 3 

enter  magnification 

• 1 

l*color,  2=erase,  3=magnify,  4=window,  6=file  print /plot,  6=transform 

• 0 

l=triangles,  2=contours,  3=surface,  4=profiles,  5=flow  lines,  6=settings 

• 1 

enter  level 

• 1 

plot  0 . 360  sec . 

enter  level 

• 2 

plot  0.620  sec. 

enter  level 

• 0 

l=triangles,  2=contours,  3=surface,  4=profiles,  5=flow  lines,  6=settings 

• 2 

pde  no.,  l=f unction/2='abs (grad) , l=fill/2=no  fill 
112 

number  of  contours 
• 11 

min  and  max  values  O.OOOOOOE+00  l.OOOOOOE+00 
contour  values 

8.333E-02  1.667E-01  2.600E-01  3.333E-01  4.167E-01 

6.000E-01  6.833E-01  6.667E-01  7.600E-01  8.333E-01 

9.167E-01 

plot  0 . 480  sec . 

pde  no.,  l=function/2=abs(grad) , l=fill/2=no  fill 

• 0 0 0 

l=triangles,  2=contours,  3=surface,  4=profiles,  5=flow  lines,  6=settings 

• 3 

pde  no. , l=function/2=abs(grad) 

• 00 

l=triangles,  2=contours,  3=8tirfac8,  4=profiles,  6=flow  lines,  6=settings 

• 5 

pde  no.,  number  of  points,  number  of  flow  lines 

• 1 2 11 

enter  2 x-values  followed  by  2 y-values 

• "•S  "•S  ,6 

total  +,  - flow  9.743972E-01  O.OOOOOOE+00 
plot  0.760  sec. 

pde  no.,  number  of  points,  number  of  flow  lines 

• 0 0 0 
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l=triangles , 2=contours,  3=surface,  4=profiles,  6=flow  lines,  6=setting3 

• 4 

pde  no.,  number  of  profiles 

• 1 21 

indices  of  profile  planes 
• 10 

plot  0.910  sec. 
pde  no.,  number  of  profiles 
• 00 

l=triangles,  2=contours,  3=surface,  4=profiles,  5=flow  lines,  6=settings 

• 0 

l=print,  2=plot,  3=iter,  4=refine,  5=comp,  6=save,  7=diff,  8=end,  9=change 

• 8 
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Appendix  D.  Program  to  Solve  Nonlinear  System  Example 

external  aex,  fex,  gex 
integer  io(9) 

parameter  (irsize  = 400000) 

parameter  (idsize  = 200000) 

c ommon/c  St  ak/dst  eik 

double  precision  dstak(idsize) 

real  rstak(irsize) 

equivalence  (dstak(l) , rstak(l)) 

save  /cstak/ 

call  istkinCirsize . 3) 

io(l)  » 5 

io(2)  = 6 

io(3)  » 125 

io(4)  * 20 

io(5)  * 21 

io(6)  » 22 

io(7)  » 23 

io(8)  = 24 

io(9)  = 26 

do  10  j * 1,  9 

call  iosetCj,  io(j)) 

10  continue 

openCunit  = io(4) , file  = ’input2.dat’,  status  = ’old*) 
openCunit  * io(5),  file  = ’plot.dat’  , status  = ’unknown’) 
openCunit  = io(6) , file  = ’print.dat’  , status  = ’unknown’) 
if  ( io(7)  .gt.  0)  then 

openCunit  = ioC7) , file  = ’insol.dat’,  status  = ’unknown’, 

♦ form* ’unformatted’) 
else 

openCunit  * -ioC7) , file  = ’insol.dat’,  status  = ’unknown’, 

♦ form* ’formatted’) 
endif 

if  C ioC8)  .gt.  0)  then 

openCunit  * ioC8) , file  = ’outsol.dat’,  status  = ’unknown’, 

♦ form* ’unformatted’) 
else 

openCunit  * -ioC8) , file  = ’outsol.dat’,  status  = ’unknown’, 

♦ form* ’formatted’) 
endif 

openCunit  * ioC9) , file  * ’scrtch.dat’,  status  * ’unknown* , 

♦ form* ’unformatted’) 
call  drvallCaex.  fex.  gex) 
closeCioC4) ) 
closeCioCS)) 
closeCioCd) ) 

close CiabsCioC7))) 
closeCiabsCioC8) ) ) 
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close (io (9) ) 

stop 

end 

subroutine  aexCnpde,  ndim,  dopde,  x,  y.  i,  u,  dudx,  dudy, 

* a.  dadu,  dadux.  daduy) 
real  x,  y,  u(*) , dudx(*) , dudy(*) , 

♦ a(*) , daduCndim,*) , daduxCndim,*) , daduy (ndim, *) 
logical  dopde (npde) 

do  10  j =»  1 , npde 
do  10  k = 1 , npde 
dadu  (j,  k)  = 0.0 
daduxCj , k)  =0.0 
daduy ( j , k)  = 0.0 
10  continue 

a(l)  = 1.0  + dudx(2)**2  + dudy(2)**2  + dudx(3)**2  + dudy(3)**2 

daduxCl,  2)  = 2.0*dudx(2) 

daduy (1,  2)  = 2.0*dudy(2) 

daduxCl,  3)  = 2.0*dudx(3) 

daduy (1.  3)  = 2.0*dudy(3) 

a(2)  = 1.0 

a(3)  » u(l)*u(2) 

dadu(3,  1)  = u(2) 

dadu(3,  2)  = u(l) 

ret\im 

end 

subroutine  f ex (npde,  ndim,  dopde,  x,  y,  i,  u,  dudx,  dudy. 

♦ f,  dfdu,  dfdiox,  dfduy) 
real  x,  y,  u(*) , dudx(*) , dudy(*) , 

* f(*),  dfdu(ndim, *) , df dux (ndim, *) , dfduy (ndim, *) 
logical  dopde (npde) 

do  10  j =1,  npde 
do  10  k = 1 . npde 
dfdu  (j,  k)  = 0.0 
dfdux(j , k)  =0.0 
dfduy (j , k)  =0.0 
10  continue 

f(l)  = 2.0*u(l) 

dfdud,  1)  * 2.0 

f(2)  = u(2)**2  * exp(-x) 

dfdu(2,  2)  = 2.0  ♦ u(2)  ♦ exp(-x) 

f(3)  = 0.0 

return 

end 

subroutine  gex(npde,  ndim,  dopde,  x,  y,  i,  j,  u,  g,  dgdu) 
real  x,  y,  u(*) , g(*) , dgdu(ndim,*) 
logical  dopde (npde) 
g(l)  = u(l)*u(2)*u(3)  - 1.0 
dgdud,  1)  = u(2)*u(3) 
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dgdu(l,  2)  = u(l)*u(3) 
dgdu(l,  3)  = u(l)>«'u(2) 
g(2)  * u(2)  - exp(x) 
dgdu(2,  2)  = 1.0 
g(3)  = u(3)  - exp(y) 
dgdu(3,  3)  = 1.0 
return 
end 

subroutine  usrgssCnpde,  x.  y,  i,  u) 
real  u(*) 
do  10  n = 1,  npde 
u(n)  =1.0 
10  continue 
return 
end 

subroutine  usrshf (ixtype , npde,  nv,  u,  xmew,  vx,  vy) 
real  u(npde,nv),  unew(nv) , vx(nv) , vy(nv) 
if  (ixtype  .le.  npde)  then 
nth  = ixtype 
do  10  i = 1,  nv 

unew(i)  = u(nth,  i)  - utrue(nth,  vx(i) , vy(i)) 
10  continue 

else 

do  20  i = 1 , nv 

unew(i)  = u(l,i)  ♦ u(2,i)  * u(3,i) 

20  continue 

endif 
return 
end 

real  function  utrueCnth,  x,  y) 
real  x,  y 

if  (nth  .eq.  1)  utrue  = exp(-x-y) 

if  (nth  .eq.  2)  utrue  » exp(x) 

if  (nth  .eq.  3)  utrue  » exp(y) 

return 
end 
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Appendix  E.  File  “input2.dat’*  for  Nonlinear  System  Example 
It  two  squares 

It  dirichlet  b.c.  on  all  sides 

0 It  old  solution  flag 

It  nt  nv  nc  nb  npde 

4 6 0 6 3 

# vertices 

# i Ivl  X y 

1 1 -1.0  1.0 
2 1 0.0  1.0 

3 1 1.0  1.0 

4 3 -1.0  0.0 

6 1 0.0  0.0 
6 1 1.0  0.0 

# curved  edge  centers 

# (none) 

# triangles 

# i vl  v2  v3 
112  4 

2 2 3 5 

3 2 4 5 

4 3 5 6 

# boundary  conditions: 

# i triangle  edge  curve  b.c. 

112  0 111 

2 13  0 111 

3 2 3 0 1 1 1 

4 4 2 0 1 1 1 

5 4 10  111 

6 3 10  111 

2 8 0 # maxm  maxr  ittype 

407  # Ivlref  Ivlusr  Ivlmax 

20  # threshold  percentage 

1 # print  level 

350000  # work  space 
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Appendix  F.  Transcript  of  interactive  session  for  Nonlinear  System  Example 
Output  from  the  computer  is  in  typewriter  font.  Numbers  typed  by  the  user  are  in 
boldface  and  have  a • at  the  beginning  of  the  line. 


refsol  entered  at  level  1 6 vertices 

grid  16  vertices 

l=print,  2=plot,  3=iter,  4=refine,  5=comp,  6=save , 7=diff,  8=end,  9=change 

• 3 


enter  niter,  ieps,  relerr,  abserr,  jsmaoc,  param 
• 10  0 0 .1  5 0 
Iter  1 


1 newton 

nom.  solution 

norm,  relerr 

3 . 437E+00 

l.OOOE+00 

3.437E+00 

2 newton 

norm. solution 

norm, relerr 

1.718E+00 

l.OOOE+00 

1.718E+00 

3 newton  norm. solution  norm, relerr 
using  newton  step  fraction  2.500E- 

1.718E+00 

•01 

l.OOOE+00 

1.718E+00 

Iter  2 

1 newton 

norm,  solution 

norm, relerr 

6.817E-01 

1 . 158E+00 

5.887E-01 

2 newton 

norm,  solution 

nom,  relerr 

1.289E+00 

1.430E+00 

9.015E-01 

3 newton 

norm. solution 

nom,  relerr 

1.289E+00 

1.430E+00 

9.015E-01 

Iter  3 

1 newton 

norm,  solution 

nom,  relerr 

8.786E-01 

1 . 840E+00 

4.775E-01 

2 newton 

norm, solution 

nom,  relerr 

8.669E-03 

2.718E+00 

3.189E-03 

3 newton 
Iter  4 

norm. solution 

nom.  relerr 

5.293E-02 

2.718E+00 

1.947E-02 

1 newton 

norm,  solution 

nom,  relerr 

6.897E-03 

2.718E+00 

2.169E-03 

2 newton 

norm,  solution 

nom,  relerr 

4.080E-06 

2.718E+00 

1.501E-06 

3 newton 

norm, solution 

nom,  relerr 

1.173E-02 

2.718E+00 

4.315E-03 

l=print,  2=plot,  3=iter,  4=refine,  5=conqp,  6=save,  7=diff,  8=end,  9=change 

• 4 


refsol  re-entered  at  level 

1 estimated  maximum  error 

2 estimated  maximum  error 

3 estimated  maximum  error 


2 16  vertices 

7.12E-02 
6.13E-02 
6.38E-02 


estimated  individual  PDE  errors 


0=continue,  l=print,  2=plot 

• 0 

estimated  total  error  (as  PDE  #1) 
0=continue , l=print , 2»plot 

• 0 


refine  to  level  2 27  vertices 

l=print,  2=plot,  3=iter,  4=refine,  5=comp,  6=save,  7=dif f , 8=end,  9=change 

• 5 


compressed  to  level  1 

l=print,  2=plot,  3=iter,  4=refine,  6=comp,  6=save,  7=dif f , 8=end,  9=change 

• 3 

enter  niter,  ieps.  relerr.  abserr,  jsmax.  param 

• 10  0 0 .01  5 0 

Iter  1 

1 newton  norm, solution  norm, relerr  2.014E-01  2.718E+00  7.411E-02 
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2 

newton 

norm, solution 

norm, relerr 

2.14©E-©1 

2.718E+00 

7.874E-©2 

3 newton 
Iter  2 

norm, solution 

norm, relerr 

2.1©4E-©1 

2.718E+©© 

7.741E-©2 

1 

newton 

norm, solution 

norm, relerr 

l.©©7E-©2 

2.718E+©© 

3.7©6E-©3 

2 

newton 

norm,  solution 

norm,  relerr 

9.477E-©4 

2.718E+©© 

3.486E-©4 

3 newton 
Iter  3 

norm,  solution 

norm, relerr 

3.569E-©2 

2.718E+©© 

1.313E-©2 

1 

newton 

norm,  solution 

norm, relerr 

7 . 878E-©6 

2.718E+©© 

2.898E-©6 

2 

newton 

norm, solution 

norm, relerr 

4.161E-©8 

2.718E+©© 

1.527E-©8 

3 

newton 

norm, solution 

norm, relerr 

1.327E-©4 

2.718E+©© 

4.883E-©5 

l=print,  2=plot,  3=iter,  4=refine,  5=comp,  6=save,  7=diff,  8=end.  9=change 

• 4 


ref sol  re-entered  at  level 

1 estimated  maximum  error 

2 estimated  maximiim  error 

3 estimated  meacimum  error 


2 27  vertices 

1.70E-02 
6.18E-02 
7.46E-02 


estimated  individual  PDE  errors 


©“continue , 1 “print , 2=plot 

• 0 


estimated  total  error  (as  PDE  #1) 
©“continue , 1 “print , 2=plot 

• 0 


refine  to  level  2 45  vertices 

Imprint,  2*plot,  3“iter,  4=refine,  6=comp,  6=save , 7»diff,  8=*end,  9*change 

• 3 

enter  niter,  ieps,  relerr,  abserr,  jsmax,  param 

• 10  10  0 .001  5 0 
Iter  1 


Chkmtx  136  equations.  Worst  equation  1©  9.894E-©2 

This  set  of  equations  may  be  unsuitable  for  Gauss-Seidel  iteration. 

1 newton  norm, solution  norm, relerr  7.868E-©2  2.718E+©©  2.891E-©2 

2 newton  norm, solution  norm, relerr  6.731E-©2  2.718E+©©  2.476E-©2 

3 nevrton  norm, solution  norm, relerr  6.66©E-©2  2.718E-*-©©  2.446E-©2 

Iter  2 

Chkmtx  135  equations.  Worst  equation  1©  9.894E-©2 

This  set  of  equations  may  be  unsuitable  for  Gauss-Seidel  iteration. 

1 newton  norm, solution  norm, relerr  1.728E-©3  2.718E+©©  6.357E-©4 

2 newton  norm, solution  norm, relerr  9.©14E-©6  2.718E-*-©©  3.316E-©6 

3 newton  norm, solution  norm, relerr  4.439E-©3  2.718E-*-©©  1.633E-©3 

Iter  3 

Chkmtx  135  equations.  Worst  equation  1©  9.894E-©2 

This  set  of  equations  may  be  unsuitable  for  Gauss-Seidel  iteration. 

1 newton  norm, solution  norm, relerr  6.482E-©6  2.718E+©©  2.385E-©6 

2 newton  norm, solution  norm, relerr  1.491E-©7  2.718E+©©  5.485E-©8 

3 newton  norm, solution  norm, relerr  1.462E-©6  2.718E-*-©©  5.343E-06 

l“print,  2“plot,  3=iter,  4“refine,  6“comp,  6“save , 7*diff,  8»end,  9-change 
• 1 

l=8\innnary , 2“timing,  3=storage,  4=vertices,  5=pde  variables,  6-settings 


52 


• 1 


n 

Ivl 

vertices 

triangles 

unorm 

enorm 

digits 

1 

1 

27 

50 

9.40E-01 

3.25E-02 

1.46 

2 

1 

27 

50 

1.40E+00 

6.81E-02 

1.31 

3 

1 

27 

50 

1.84E+00 

7.25E-02 

1.40 

1 

2 

45 

90 

2 

2 

45 

90 

3 

2 

45 

90 

l=summary,  2=timing,  3=storage,  4=verbices,  6=pde  variables,  6=settings 

• 2 

timing  distribution 
seconds  percent 


plot  0.000  0.00 
assembly  7.980  50.44 
factor  2.480  15.68 
mg  1.210  7.65 
rhs  check  3.080  19.47 
refine  1.070  6.76 
total  15.820  100.00 


l=8ummary . 2*timing,  3=storage,  4=vertices,  5=pde  variables,  6=settings 

• 3 

45  vertices  used. 

1710  vertices  allowed. 

97  trieingles  used. 

4559  triangles  allowed. 

l=summary , 2-timing,  3=storage,  4=vertices,  5=pde  variables,  6=settings 

• 0 

l=print,  2=plot,  3=iter,  4=refine,  5=corap,  6=save,  7=dif f , 8=end,  9=chaoige 

• 2 

l^triamgles,  2*contours,  3=surface,  4=profiles,  5=flow  lines,  6=settings 

• 1 

enter  level 

• 1 

plot  0.150  sec. 

enter  level 

• 2 

plot  0 . 220  sec . 

enter  level 

• 0 

l=triauigles , 2=contours,  3=surface,  4*profiles,  5=flow  lines,  6=settings 

• 6 

l=color,  2=erase,  3=magnify,  4=window,  5=file  print/plot,  6=transform 

• 6 

enter  index  of  transformation 

• 2 

transformation  number  2 as  PDE  number  1 

l=color,  2=erase,  3=magnify,  4=window,  5=file  print /plot,  6=transform 

• 0 
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l=triangles,  2=contours,  3=surface,  4=profiles,  5=flow  lines,  6=settings 

• 3 

pde  no.,  l=f\mction/2=abs(grad) 

• 11 

level,  l=surface  triangles,  2=surface  contours,  3=surface  grids 
• 21 

viewing  direction  — nx,  ny,  nz 
• 1-21 

min  and  max  values  -1 .334548E-03  O.OOOOOOE+00 

plot  0.940  sec. 

pde  no.,  l=function/2=abs(grad) 

• 11 

level,  l=surface  triangles,  2=surface  contours,  3=surface  grids 

• -2  2 

number  of  contours 

• 9 

viewing  direction  — nx,  ny,  nz 

• 1-21 

min  and  max  values  -1 .334648E-03  O.OOOOOOE+00 
contour  values 

-1.201E-03  -1.068E-03  -9.342E-04  -8.007E-04  -6.673E-04 
-5.338E-04  -4.004E-04  -2.669E-04  -1.336E-04 
plot  0.460  sec. 

pde  no.,  l=function/2=abs(grad) 

• 11 

level,  l=surface  triangles,  2=surface  conto\irs,  3=surface  grids 

• 23 

number  of  x lines,  number  of  y lines 

• 21  11 

viewing  direction  — nx,  ny,  nz 
• 1-21 

min  and  max  values  -1 .334648E-03  O.OOOOOOE+00 

plot  1.610  sec. 

pde  no.,  l*function/2=abs(grad) 

• 00 

l=triangles,  2=contours,  3=surface,  4*profiles,  5=flow  lines,  6=settings 

• 0 

Imprint,  2=plot,  3=iter,  4=refine,  6=comp,  6=8ave,  7=diff,  8=end,  9»change 

• 8 
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Figure  3 


Figure  4 
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